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CHAPTER  1 


ENERGY  FLOW  IN  A  DEGENERATE  SEMICONDUCTOR 
DEVICE  WITH  NONUNIFORM  BAND  STRUCTURE 

AND  RADIATION 


1.1  Introduction 

The  simulation  of  heat  or  energy  flow  in  semiconductor  devices  is  a  complex  problem. 
There  are  three  basic  models  for  energy  transport:  The  first  one  is  due  to  Stratton  [1], 
which  represents  a  good  approximation  as  long  as  the  asymmetric  portion  of  the  momen¬ 
tum  distribution  of  carriers  is  small.  The  so-called  hydrodynamic  model  due  to  Bfotekjaer 
[2]  is  slightly  more  general  but  has  the  disadvantage  that  several  transport  coefficients 
can  only  be  determined  by  direct  solution  of  the  Boltzmann  transport  equations  (BTE). 
The  third  is,  of  course,  the  direct  solution  of  the  BTE  by  the  Monte-Carlo  methods,  for 
example.  Here  we  use  Stratton’s  approach  with  a  careful  choice  of  proper  expressions  for 
particle  and  energy  or  heat  fluxes.  Expressions  for  particle  fluxes  derived  from  the  BTE 
have  been  studied  extensively  [3]— [12]. 

We  begin  with  a  discussion  of  the  coupling  of  the  energy  conservation  equations 
for  various  physical  subsystems,  including  the  radiation  system,  as  appropriate  for  the 
simulation  of  optoelectronic  semiconductor  devices.  We  examine  various  mathemati¬ 
cally  equivalent  expressions  for  their  numerical  performance  in  actual  simulations.  The 
discretized  heat  equation  often  shows  signs  of  an  ill-conditioned  matrix  problem.  We 
investigate  causes  and  cures,  and  discuss  the  discretization  scheme  for  various  flux  ex¬ 
pressions. 

We  consider  three  types  of  boundary  conditions  for  the  thermal  equation:  Dirichlet 
conditions,  von  Neumann  conditions  with  vanishing  energy  flux,  and  heat  current  bound¬ 
ary  condition.  In  a  realistic  situation,  lead  wires  are  typically  used  to  provide  a  current 
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or  voltage  bias.  We  also  review  the  adaptability  of  the  various  equivalent  expressions  for 
particle  and  energy  fluxes,  and  attempt  to  clarify  physical  concepts  of  heat  flux  and  dis¬ 
sipative  energy  flux  in  the  energy  transport  equation.  These  concepts  are  often  borrowed 
from  classical  nonequilibrium  thermodynamics  and  need  to  be  adapted  to  the  electron 
dynamics  in  semiconductors. 

1.2  Conservation  of  Energy 

We  consider  four  physical  systems:  the  systems  of  electrons,  the  holes,  the  crystal 
lattice,  and  the  radiation.  We  then  can  construct  four  energy  conservation  equations  for 


these  four  systems  (Appendix  A): 

dujdt  =  —V  •  Se  —  EqUhSR  ~  Qe-lat 

+  ^Aug  +  -^e-h  ~  ^C,tuZ^ rad  +  huR^t fc,  (1-1) 

du^/dt  =  -V  •  Sh  -  E^Ursr  -  Qh-iat 

+  AAug  4-  Ae-h  —  Ey  h-UTad  +  ha>f?h,fc)  (1-2) 

duujdt  =  -V  -  Su*  +  EfUns  r  +  Q,  (1.3) 

durad/dt  =  -V  .  Srad  +  hu  (i UTad  -  Rlc)  ,  (1.4) 

E+  —  Eq  —  Ev ,  hu  =  Ec'Kz  —  EVK3 , 

Rtc  =  R{:  -  Rt  Q  =  ge-lat  -  Qh-lat. 


Here  t  is  time,  and  tte,  Uh,  uut,  and  uraa  are  the  energy  densities  for  the  electron,  the 
hole,  the  lattice,  and  the  radiation  systems,  respectively.  The  vectors,  Se,  Sh,  Siat, 
and  Srad  are  the  energy  flux  densities  for  the  respective  systems.  The  energy  flux  for 
the  radiation  system,  Sraa,  is  also  called  the  Poynting  vector  in  electromagnetics.  The 
symbols,  UT& a  and  C/hsr,  represent  the  net  rates  of  radiative  recombination  and  Hall- 
Shockley-Reed  recombination,  and  R *c  and  /?h  refer  to  the  free-carrier-absorption  rate 
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for  electrons  and  holes,  respectively.  The  energy  denoted  by  E+ ,  which  is  slightly  larger 
than  the  energy  band  gap  Es,  is  the  average  energy  released  by  the  Hall-Shockley-Reed 
recombination  or  the  Auger  recombination  depending  on  where  it  is  used.  The  constant, 
h  =  h/2ir,  is  the  Planck  constant,  and  hU  is  the  average  energy  released  by  the  radiative 
recombination.  The  terms  XaU8  and  Ae-h  represent  the  rate  of  energy  exchange  due 
to  the  Auger  recombination  and  electron-hole  scattering,  respectively.  The  terms  Qc 
and  Qh  represent  the  heat  dissipation  to  the  lattice  system  from  the  respective  carrier 
systems  by  the  in.raband  thermal  relaxation  process.  The  sources  of  these  terms  are  the 
Joule  heat,  the  Peltier  heat  (at  abrupt  heterojunctions),  and  the  heat  dissipated  from 
excited  electrons  or  holes  due  to  free  carrier  absorption  and  the  Auger  process.  The 
Auger  process  itself  does  not  involve  energy  transfer  to  the  lattice  system.  However,  the 
Auger  recombination  process,  a s  well  as  the  free-carrier  absorption  process,  leaves  an 
electron  in  mid-conduction  band  or  a  hole  in  mid- valence  band.  This  carrier  then  relaxes 
its  energy,  giving  off  heat  to  the  lattice.  The  energy  transferred  to  the  lattice  system  by 
this  process,  as  well  as  the  Joule  heat,  is  included  in  the  term  Q. 

The  radiation  system  and  the  lattice  system  can  be  thought  of  as  sums  of  energy 
conservation  equations  for  infinitely  many  modes.  The  radiation  system  is  often  conve¬ 
niently  represented  by  the  photon  rate  equations  for  individual  modes  after  integrating 
over  the  space  under  consideration.  Each  photon  rate  equation  is  of  the  form 

dSJdt  =  G„SV  +  R?  -  5„/rph.  (1.5) 

Here  t„  is  the  photon  lifetime,  and  Sv,  Gu  and  are  the  photon  occupation  number,  the 
mode  gain,  and  the  spontaneous  emission  rate  of  mode  v.  Note  that  the  space  integration 
of  the  Poynting  vector  has,  with  that  of  the  free-carrier  absorption  term,  been  replaced 
by  S„/tu,  and  that  of  the  optical  energy  density  has  been  replaced  by  the  sum  of  S„'s. 

The  lattice  system  cannot  be  treated  in  this  way.  The  anharmonicity  of  the  Hamilto¬ 
nian  creates  a  coupling  of  phonon  modes  resulting  in  thermodynamically  irreversible  heat 
dissipation  through  Umklapp  processes.  Energy  exchange  between  the  carrier  system  and 
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the  resonant  radiation  system  (as  in  a  semiconductor  laser)  maintains  the  reversibility 
on  short-time  scales,  resulting  in  transient  “ringing”  of  radiation  mode  amplitudes. 

The  information  which  we  seek  to  obtain  from  the  heat  flow  simulation  is  the  in¬ 
ternal  temperature  distribution.  If  the  energy  exchange  among  the  various  systems  is 
weak,  one  is  forced  to  introduce  individual  average  energies  for  electrons,  holes,  lat¬ 
tices,  and  photons.  If  the  energy  exchange  between  the  electron  system  and  the  hole 
system  is  fast  enough,  the  average  electron  and  hole  energies  can  be  considered  equal, 
and  we  need  only  one  conservation  equation  for  those  two  systems  with  the  source  term 
-Ku  (Utad  —  i?fc)  —  E+(JHsn  —  Q.  If  the  energy  exchange  between  the  carrier  and 
lattice  systems  is  rapid,  we  can  introduce  a  temperature  which  is  equal  for  electrons, 
holes,  and  the  crystal  lattice,  and  we  may  combine  (1.1)— (1.3)  to  obtain 

£  (ue  -  Uh  +  «u0  =  -V  ■  (Se  -  Sh  +  Slat)  -  fiw  (f/rad  -  Rk)  .  (1.6) 

In  doing  this  we  have  canceled  the  heat  dissipation  term  Q  for  the  Joule  heat  and  the 
hot  carrier  relaxation,  either  of  which  is  not  easy  to  realize  or  have  no  simple  explicit 
expression  for  numerical  simulation. 

In  the  simulations  we  also  need  to  include  the  Poisson  equation  and  the  following 
carrier  continuity  equations  along  with  energy  balance  equations. 

Q 

—  {n,  p}  +  V  •  {je,  jh}  +  UTad  +  C^HSR  +  t^Aug  =  0.  (1.7) 

Here  n  [p]  and  je  [jh]  are  the  electron  [hole]  density  and  the  particle  flux  density,  re¬ 
spectively,  and  f/Aug  is  the  net  rate  of  the  Auger  recombination.  To  solve  all  equations 
self-consistently,  one  needs  to  express  the  various  contributions  to  those  equations  in 
terms  of  a  set  of  variables  which  are  chosen  as  the  unknowns  for  the  system  of  equations. 
These  variables  are  the  electrostatic  potential,  the  carrier  densities  n  and  p,  and  the 
temperature,  under  the  assumption  that  local  equilibrium  is  maintained  so  that  a  local 
temperature  is  a  valid  concept.  Instead  of  n  and  p,  we  use  r/e  and  rj^,  sometimes  called 
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“Planck  potentials”  [13]  which  are  defined  below.  This  is  advantageous  when  the  heat 
transfer  is  involved. 


1.3  Analytic  Considerations 

1.3.1  Expressions  for  the  conserved  quantities  and  sources 

There  are  various  expressions  available  for  the  various  rate  terms  in  (1.5)— (1.7)  with 
varying  degrees  of  approximation,  which  are  described,  e.g.,  in  [14]. 

The  conserved  quantities  in  the  continuity  equations  can  be  expressed  in  terms  of  rje, 
%,  and  temperature.  Assuming  a  parabolic  band  structure  with  density-of-states-equi- 
valent  effective  masses  me  and  for  electrons  and  holes,  respectively,  one  has 


«  =  NkT^in(V.),  p  = 

(1.8) 

*  iV£T3/I  [jT*,/^.)  +  EcFy /,(>,.)] , 

(1.9) 

u„  ~  KT1'1  [jT-F ,/,(%)  -  ■Ev.F.M-fc)] , 

(1.10) 

/q  =  2  (me/2jrfi2)  ^  ,  JVy  =  2  (mi,/2Th2j  *  , 

p,  =  (F,-Ec)/T,  ^  =  (Ey  -  Fb) /T, 

(1.11) 

£ 

II 

1 

►O 

1 

>< 

S3 

II 

£ 

1 

o 

(1.12) 

l  r00  x1 
o  e*-*  +  1 


dr, 


where  Fe  and  F b  are  the  quasi-Fermi  levels  for  electrons  and  holes,  respectively,  Eq  and 
Ev  are  the  conduction  and  valence  band  edges,  respectively,  and  xp  is  the  electrostatic 
potential.  The  variables  x  and  Eq  refer  to  the  electron  affinity  and  bandgap  of  the 
material,  respectively,  and  T  is  the  temperature  in  energy  units.  Closed-form  expressions 
for  the  above  quantities  for  nonparabolic  band  structure  axe  found  in  [7].  The  effect  of 
particle  streaming  on  the  carrier  densities  has  been  included  within  the  framework  of  local 
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quasi-Fermi  levels  and  the  local  temperatures.  However,  the  effect  on  the  energy  densities 
2melie|2/n5  which  is  additional  to  (1.9),  has  been  ignored  by  taking  the  equilibrium 
Fermi-Dirac  functions  for  the  integration  over  the  phase  space.  This  should  be  of  concern 
as  discussed  in  Appendix  A  regarding  the  equation  for  the  conservation  of  energy. 

The  energy  stored  in  the  lattice  is  [9] 


uut 


D{») 

eh*/T  _  i 


do;, 


(1.13) 


where  u?d  represents  the  Debye  frequency,  and  D(ut)  is  the  density  of  phonon  states.  For 
use  in  the  context  of  time  derivative,  the  following  equation  is  used  [13]. 


duxtx/dt  =  cpdTjdt, 


(1.14) 


where  cp  is  the  pressure-specific  heat  capacity  per  unit  volume  of  the  device  material. 

The  expression  for  each  term  of  hu  (Urad  —  R{c^j  in  (1.6),  which  represents  the  energy 
exchange  from  the  stimulated  emission,  spontaneous  emission,  and  free-carrier  absorp¬ 
tion,  are  discussed  in  Chapter  2. 


1.3.2  Expressions  for  fluxes 

In  addition  to  the  above  rates,  expressions  for  the  various  fluxes  are  needed.  We 
assume  that  the  relaxation  time  of  hot  carriers  is  negligibly  small  on  our  time  scale. 
Then  we  can  use  the  linear  form  of  the  transport  equations  derived  from  the  BTE  for 
the  particle  fluxes  in  the  absence  of  a  magnetic  field  in  a  degenerate  semiconductor.  The 
particular  form  of  expressions  for  various  fluxes  which  we  use  here  is  most  suitable  for 
the  case  when  the  carrier  degeneracy,  band  structure  nonuniformity,  and  temperature 
gradient  are  important  at  the  same  time. 

je  =  — M„  [jF0(T?e)  (TVr,e  +  VEC)  +  2^1(7?e)VT] ,  (1.15) 

Jh  =  -A4  [^ofo.)  (TV;*,  -  VEv)  +  2^,(7 7h)VT] ,  (1.16) 
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Me  =  /zen/g^0(j?e),  .Mh  =  fihp/q? r0(»7h).  (1-17) 

Here  pe  and  ph  are  the  mobilities  of  electrons  and  holes,  respectively,  and  <7  is  the  el¬ 
ementary  charge.  The  energy  flux  expressions  which  are  symmetric  in  the  sense  of  the 
Onsager  relations  [13]  are 

s.  =  Eci.  +  S'r,  Sh  =  Bvik  -  sf1,  (1.18) 

sf  =  -M.r(2^i(i).)  (TV,.  +  V£c)  +  6^(r,.)VT] ,  (1.19) 

Sh“  =  (TV%  -  VEv)  +  6F2(v,)VT]  .  (1.20) 

The  derivation  of  the  above  particular  forms  with  an  extension  to  the  first-order  approxi¬ 
mation  for  the  nonparabolic  bands  is  given  in  Appendix  B,  along  with  integral  expressions 
deduced  from  the  BTE. 

The  heat  flux  equation  for  phonons, 


Sut  =  (1.21) 

can  also  be  used  as  the  expression  for  the  phonon  energy  flux,  with  k  being  the  lattice 
heat  conductivity  of  the  material. 


1.3.3  Quantum-well  regions 

The  expressions  in  (1.8)— (1 . 12)  are  used  everywhere  except  in  the  size-quantized  re¬ 
gions.  Inside  the  quantum  well,  we  use 
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where  and  axe  the  heavy  and  light  hole  masses  in  the  quantum  well,  respectively, 
which  may  differ  from  the  values  for  a  bulk  semiconductor.  In  this  way,  we  have  been 
able  to  simplify  expressions  for  various  physical  quantities  when  the  temperature  is  also 
a  variable  in  the  system  of  equations. 

The  particle  current  at  the  abrupt  hetero-interface  is  due  to  thermionic  emission.  The 
net  electron  current  toward  the  outside  of  the  quantum  well  at  the  interface  is  simulated 
as 


U  = 


—  ■ -rpl  r 
Tfl+J 


2 


/in  SC*  -  E'c\  t 

expU* - - J-exp(T/e  j 


m=2/  (l/m“  +  l/meout)  ,  T  =  (fin  +  Tout)  /2, 


where  is  the  effective  mass  of  electrons  for  the  node  at  the  quantum-well  side  (lower 
Ec),  while  m°ut  is  that  for  the  side  of  higher  band  edge  Ec .  Variables  such  as  j/J",  r?°ut, 
Eq  ,  and  Eq11  are  all  similarly  defined.  Also,  the  corresponding  energy  current  should  be 
simulated  accordingly.  That  is, 


Se  =  (££ut  4-  Tm  +  Tout)  je 


at  the  hetero- interface.  The  net  hole  current  and  hole  energy  current  are  similarly  ex¬ 
pressed,  taking  both  the  heavy  and  light  holes  into  account.  Previous  simulations  of  het¬ 
erostructure  lasers  [15]— [19]  were  based  solely  on  the  drift-diffusion  model  for  the  carrier 
transport.  This  model  yielding  over-estimated  drift  velocity  at  the  abrupt  hetero-inter¬ 
face  is  not  only  physically  incorrect,  but  also  prone  to  poor  convergence  if  the  electrostatic 
potential  varies  strongly  around  the  abrupt  hetero- interface  with  high  injection  currents. 


1.4  Conditioning  the  Heat  Flow  Equation 

In  numerical  simulations  with  finite  floating-point  precision,  one  needs  to  consider 
the  degree  of  conditioning  of  the  system  of  partial  differential  equations  (in  our  case, 
(1.6)— (1.7)).  Notice  that  je  and  jh  are  multiplied  by  Ec  and  E\  in  (1.19)  and  (1.20), 
respectively.  The  energies  Ec  and  Ey  have  been  cited  with  no  reference  level,  since  the 
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zero  point  of  energy  can  be  chosen  freely  as  long  as  it  is  consistent  throughout.  Then 
Ec  and  Ey  are  determined  by  solving  the  Poisson  equation,  and  the  absolute  values  of 
these  two  quantities  from  (1.12)  are  usually  much  greater  than  T ,  the  thermal  energy. 
Therefore,  the  magnitude  of  the  first  terms  in  (1.19)— (1.20)  is  much  greater  than  the 
remainder  of  the  equations.  As  a  consequence,  the  heat  flow  equation,  (1.6),  becomes 
approximately  equal  to  a  linear  combination  of  the  carrier  continuity  equations,  (1.7), 
whenever  the  recombination  terms  or  the  lattice  heat  flux  term  is  much  smaller  than 
\Ecje  —  Eyj h|  at  any  portion  of  the  device  profile.  This  results  in  an  ill-conditioned 
matrix  problem  after  discretization  on  a  finite-precision  digital  computer,  and  may  be 
the  main  cause  of  the  poor  convergence  in  many  heat  flow  simulation  problems  which 
do  not  include  the  second  terms  in  (1.19)— (1.20),  or  are  not  properly  conditioned.  It  is, 
therefore,  of  great  importance  to  choose  an  appropriate  energy  reference  for  Ec  and  Ey. 
A  reasonable  strategy  for  this  choice  is  to  separate  reference  levels  for  the  two  band  edges. 
The  radiative  recombination  rate  term  cancels  out  by  choosing  E™{  —  (Ec  +  Ey) /2-\-hu /2 
and  Ey{  =  (Ec  +  Ey)/ 2  —  hu/2  for  the  two  reference  levels  E™*  and  Ey{,  where  Ec  and 
Ey  are  the  average  values  of  conduction  and  valence  band  edges,  respectively.  Then, 


°v 


dT 

dt 


§ TFfM  +  (i?c  -  E°  ~  y) 


-  Kri  [ +  (£V  -  e.  +  to/2)  JjOm,)]  } 

-  V  .  kVT  +  V  •  [(£c  -E0-  hu/2)  je  +  Sf1 
+  (-Ey  +  E0-hu/2)jh-S1r] 

-  hu  (UHSR  +  £/Aus  +  Rk )  =  0,  (1.24) 

E0=(eE  +  &;)/2, 
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where  E0  is  a  constant  energy  offset  chosen  as  the  average  of  Ec  and  Ey  at  the  electrode 
nodes  from  the  boundary  condition.  This  form  has  been  found  to  work  well  for  strongly 
forward-biased  devices  with  high  carrier  and  doping  densities. 

We  think  that  the  above  method  will  work  for  most  cases.  For  high  voltage  devices, 
it  may  be  necessary  to  solve  the  equations  first  without  the  energy  flow  equation  and  get 
the  profiles  of  Ec  and  Ey,  and  then  set  up  energy  references  which  are  all  different  for 
all  discretized  equations. 

1.5  Discretization  of  Flux  Equations 

To  discretize  the  system  of  partial  differential  equations  arising  from  (1.6)— (1-7),  we 
use  the  finite-difference  method  with  the  box  discretization  scheme  [20]  applied  to  rectan¬ 
gular  meshes.  A  well-known  method  for  this  type  of  discretization  problem  with  uniform 
temperature  was  first  introduced  by  Scharfetter  and  Gummel  [21].  This  method  has  been 
extended  for  the  nonisothermal  cases  with  simplified  expressions  for  the  nonisothermal 
currents  in  [22]— [23].  Here  we  show  the  results  for  our  expressions  in  (1 .15)— (1.20)  for 
the  fluxes. 

We  first  transform  the  expressions  for  the  flux  densities  for  electrons,  (1.15)  and  (1.19) 
into 

je  =  -Me  {Tl0(V')Vr0(v*)  +  Mk)  [V£c  +  27l(7e)VT]}  ,  (1.25) 

S'T  =  -2 MeT  {T7 i(fc)VJi(*) 

+  Ti(ife)  [V£c  +  372(7e)Vr]}  ,  (1.26) 

7 Av)  =  M- 

Note  that  all  orders  of  7j(7)  become  1  for  a  nondegenerate  carrier  system.  We  assume 
here  that  the  potential  and  temperature  are  almost  linear  in  the  inter 'a!  between  two 
nodes,  i  and  i  +  l,  located  at  positions  x  =  x,  and  x,+i,  respectively.  Then  the  discretized 
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form  for  the  electron  flux  in  the  +x  direction  is  found  to  be 


Xi 


sH  ,  =2rl(t^i!i-|>(4t.,WW) 

* '  2  2  *t+l 

—  B(— A^2j  02)^l(l7e|i+l)]  > 


(1.27) 


(1.28) 


B(A*i,0i)  =  A^/[exp(A^j/0J)  -  1], 

A^j  =  Ea+i  -  Ea  +  (i  +  1)  7j(fc)  (^+1  —  Ti) , 

=  7i-i(^)(7;+i  +  ri)/2, 

*?e  =  (^?e|i  d"  ^/e|«+l )  /2. 

Here,  a  function  B(- )  has  been  introduced  by  modifying  the  Bernoulli  function  B(£)  = 
£/  (e*  —  l)  to  accommodate  the  temperature  gradient  and  Fermi-Dirac  statistics.  It  is 
possible  to  use  Tq(-)  instead  of  ^x(-)  only  if  we  replace  B(A^2,  ©2)  with  ji(r^)B(A^2,  ©i)- 
That  is, 

=  2r(+.  -^±i-7.M[6(A*J,8,)^o(>)«li) 

'*+?  3  Xi+l  —  Xi 

-B(-A»j,e1)Jtb(«hli+i)].  (1-29) 


This  scheme  is  used  in  actual  implementation  of  the  algorithm. 

For  the  other  three  directions,  and  also  for  the  four  components  of  hole  flux  and 
energy  flux,  one  may  obtain  similar  expressions.  Note  that  a  nonuniform  band  structure 
such  as  graded  hetero junctions  in  compound  semiconductor  devices  are  all  included  in 
the  basic  variables,  without  introducing  additional  parameters,  which  was  suggested  in 
[24],  for  numerical  simulation. 
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Inside  the  quantum  well  with  the  staircase-like  density  of  states,  although  the  basic 
assumption  used  to  derive  the  linearized  carrier  transport  equations  in  (1.15)— (1.20)  is  not 
valid,  we  assume  that  those  equations  are  still  approximately  valid  with  the  adaptation 
that 

fo(r]e)  and  fo(Vp)  are  replaced  by  n/Nc  and  p/iVv,  (1.30) 

respectively,  where  {n,  p)  are  those  of  ( 1 .22)— ( 1 .23).  In  a  quantum- well  laser,  the  overall 
carrier  transport  is  hardly  affected  by  the  accuracy  of  the  transport  model  for  the  quan¬ 
tum-well  regions. 


1.6  Boundary  Conditions  for  the  Heat  Flow  Equation 

Boundary  conditions  represent  a  more  complex  problem  for  the  heat  transfer  equation 
than  for  the  ordinary  carrier  transport  equations.  Here  we  consider  only  special  cases. 
The  simplest  case  is  a  boundary  contacting  an  ideal  heat  sink.  This  results  in  a  Dirichlet 
boundary  condition.  The  next  simplest  case  is  the  von  Neumann  condition  with  vanishing 
heat  flux  for  the  air-semiconductor  boundaries.  For  this  condition,  one  neglects  the  small 
amount  of  heat  dissipated  to  the  air.  The  realistic  boundary  condition  at  the  electrode 
is  mor  complicated.  Suppose  that  the  total  amount  of  heat-flow  out  of  an  electrode 
is  given.  This  then  becomes  a  heat  current  boundary  condition  similar  to  the  carrier 
current  boundary  condition  for  the  carrier  continuity  equations. 

A  simple  schematic  band  diagram  showing  the  two  quasi-Fermi  levels  for  a  semicon¬ 
ductor  laser  diode  with  a  strong  forward-bias  is  shown  in  Figure  1.1.  Although  most 
injected  carriers  recombine  at  the  junction  region,  still  many  not-yet-recombined  carriers 
reach  the  electrodes.  This  is  represented  by  the  two  split  quasi- Fermi  levels  throughout 
the  device.  At  the  anode,  such  electrons  lose  energy  to  the  lattice  heating  the  contact 
region,  while  holes  gain  energy  cooling  the  same  region. 

One-dimensional  analysis  for  the  heat  generated  at  the  metal-semiconductor  Shottky 
contact  was  done  by  Stratton  [25].  Here  we  should  consider  both  electrons  and  holes 
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Figure  1.1:  Schematic  energy-band  diagram  for  the  quasi-Fermi  levels  in  a  forward- 
biased,  one-dimensional,  semiconductor  laser  diode. 

together,  and  provide  a  two-dimensional  formula  permitting  a  variation  of  current  along 
the  contact.  We  assume  that  the  contact  is  an  ideal  ohmic  contact.  The  electrons  or  holes 
then  meet  either  a  tunneling  barrier  or  a  band  discontinuity,  and  the  incoming  energy 
flux  to  the  contact  is  described  by 

(Se-Sh  +  Sut  +  Srad).n,  (1.31) 

where  n  is  the  outward  normal  unit  vector  at  the  contact.  We  only  have  a  sea  of  electrons 
which  are  supported  by  an  almost  uniform  Fermi  level  in  the  metal,  Fmet.  Then  the 
energy  flux  into  the  contact  area  is 

(Sr*  +  «£*).(-»).  (1.32) 

From  the  conservation  of  energy,  the  sum  of  (1.31)  and  (1.32)  should  be  zero.  Following 
[9],  the  energy  flux  carried  by  the  electrons  in  the  metal  is  Fmet  +  TJ^et,  where  J™et  is 
the  entropy  flux  in  the  metal,  and  TJ^et  is  sometimes  called  the  “dissipative  energy  flux 
[26].”  This  last  flux  is  again  [13] 

TJ£et  =  7rmetj 


^met  i 


jvr, 


(1.33) 
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where  nmet  is  the  Peltier  potential  [8].  The  coefficients  K™et  and  «met  represent  the  elec¬ 
tronic  and  lattice  thermal  conductivities,  respectively,  of  the  metal  used  for  the  electrode. 
The  dissipative  energy  flow  from  the  contact  to  the  bulk  metal  region  is 

Jkt  TJ$a  ■  nds  =  /«, 

where  integration  is  done  over  the  area  (in  the  metal  side)  of  semiconductor-metal  in¬ 
terface  denoted  by  ‘A*.’  For  the  heat  current  boundary  condition,  we  are  to  specify 
this  quantity  or  provide  necessary  relations  which  specify  this  quantity  at  last.  Then  an 
equivalent  expression  provides  the  necessary  boundary  condition: 

/  (Se  -  5h  +  Su t  +  Srad)  •  nds  -  /  Fmttj  ■  hds  =  /„,  (1.34) 

Ja~  Jaj 

where  fA-  represents  the  area  integration  in  the  semiconductor  side. 

Mote  that  an  electron  loses  a  considerable  amount  of  energy  when  equilibrating  at 
the  contact,  while  a  hole  gains  a  slight  amount.  However,  when  the  injection  level  is 
high,  the  number  of  holes  injected  at  the  anode  can  be  so  much  larger  than  that  of  the 
electrons  reaching  the  anode  that  the  overall  effect  cannot  be  determined  by  calculating 
only  the  energy  loss  of  electrons.  Similar  arguments  can  be  applied  to  the  cathode  side. 
This  time,  electrons  are  cooling  the  contact  region  as  a  result. 

By  assuming  an  ohmic  contact,  the  Fermi  level  of  the  metal  can  be  seen  as  the 
equilibrium  Fermi  level  of  the  semiconductor  just  inside  the  contact.  Then  the  heat 
current  boundary  condition  can  be  constructed  as 

-ni  h 

+  Skin  +  5ki„  +  S-m  +  Srad]  .  fids  =  jQ  (i.35) 

We  should  express  all  the  fluxes  in  the  left-hand  side  in  terms  of  chosen  independent 
variables.  With  this  expression  for  the  heat  current  boundary  condition,  we  can  assign 


/  [(E£m  -  F«*)je  -  (E? 

J  Ai 
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a  specific  amount  of  heat  flux  for  Iq ,  or  relate  this  expression  with  the  formula  which 
is  to  be  obtained  from  physical  consideration  of  the  environment  discussed  in  the  next 
chapter. 

Still  not  considered  is  the  light  absorption,  denoted  by  Srad,  at  the  electrode  or  the 
semiconductor-metal  interface.  In  direct  band-gap  semiconductors  with  excess  carriers, 
this  contribution  may  not  be  ignored.  For  the  stimulated  emission  in  a  semiconductor 
laser,  the  laser  light  does  not  hit  the  electrode  metal  directly.  Although  the  metallic 
surface  has  a  high  reflectivity  above  0.9,  the  tail  of  the  Fabry-Perot  mode  pattern  still 
intrudes  into  the  metal  to  the  skin  depth,  and  the  mode  intensity  is  so  great  when  lasing 
that  we  need  to  introduce  a  phenomenological  radiation  heat-coupling  coefficient  for  each 
electrode.  This  absorption  of  light  power  is  a  part  of  the  so-called  waveguide  loss,  and 
thus  the  value  of  the  coupling  coefficient  cannot  be  larger  than  what  the  waveguide 
loss  coefficient  implies.  For  the  portions  from  the  spontaneous  emission  and  from  the 
stimulated  emission,  we  thus  introduce  respective  power  coupling  coefficients  X;P  and  x*1 
as 

fSni-M.m!  (x^tC,  +  dV  (1.36) 

•M,  J 

The  right-hand  side  integration  is  carried  out  over  the  whole  device  volume.  The  values 
of  these  coefficients  are  considered  to  be  strongly  dependent  on  geometry. 
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CHAPTER  2 


MINILASE  AND  TWO-DIMENSIONAL  ANALYSIS  OF 
QUANTUM- WELL  LASERS 

2.1  Introduction 

Semiconductor  lasers  have  mainly  been  analyzed  by  use  of  analytical  or  simple  nu¬ 
merical  models.  Such  models,  however,  cannot  fully  account  for  the  multidimensional 
nature  of  light  and  carrier  interactions  in  various  laser  structures.  Even  for  a  simple 
wide-stripe  laser,  the  effects  of  laser  switching,  multimode  behavior,  and  particularly 
the  semiconductor  quantum-well  geometry  can  only  be  understood  by  elaborate  numer¬ 
ical  techniques.  Two-dimensional  simulators  have  previously  been  developed  by  other 
groups  [15]— [19].  Of  particular  interest  is  the  pioneering  work  of  Wilt  and  Yariv  [15]. 
Our  simulator,  “MINILASE,”  represents  an  improvement  for  the  design  of  quantum- well 
lasers  and  high-power  semiconductor  lasers.  It  can  account  for  heat  conduction,  temper¬ 
ature  distribution,  and  multiple  spectral  mode  behavior  as  well  as  other  two-dimensional 
characteristics  of  carrier  flow  and  optical  intensity  profile. 

The  design  goal  of  MINILASE  is  a  correct  implementation  of  semiconductor  physics. 
Since  it  solves  for  the  internal  temperature  distribution,  it  should  be  useful  in  long- 
wavelength  semiconductor  laser  applications,  where  the  threshold  current  exhibits  a 
strong  temperature  dependence,  as  well  as  for  power  semiconductor  laser  applications. 
Additionally,  the  fact  that  MINILASE  correctly  simulates  the  switching  response  of  the 
laser  will  be  useful  in  applications  such  as  short-range  optical  communications  in  optical 
integrated  circuits  as  well  as  high-speed  long-distance  data  communication. 

MINILASE  has  been  designed  as  follows:  For  the  electronic  part,  it  discretizes,  by  the 
finite-box-integration  method  [20],  four  partial  differential  equations  including  the  heat 
flow  equation.  It  solves  the  resulting  nonlinear  simultaneous  equations  by  the  full-Newton 
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method  for  both  steady-state  and  transient  responses.  For  the  optical  part,  it  solves 
the  Helmholtz  eigenvalue  equation  on  a  two-dimensional  cross-sectional  plane,  so  that 
MINILASE  can  accommodate  various  profiles  of  compound  semiconductor  compositions. 
From  an  outer  iteration  with  the  photon  rate  equation,  MINILASE  yields  the  individual 
photon  numbers  of  multiple  spectral  modes,  which  permits  spectral  mode  analysis  of  the 
optical  transition  rates  in  semiconductor  lasers. 

We  have  applied  our  simulator  to  GaAs-AlGaAs  graded-index-separate-confinemem- 
heterostructure  buried-quantum-well  (GRINSCH  buried-QW)  lasers  [27].  The  results  are 
in  good  agreement  with  the  available  experimental  data. 

2.2  Physical  Considerations 

At  least  four  equations  are  necessary  to  simulate  semiconductor  lasers;  the  Poisson 
equation  for  the  potential,  the  two  continuity  equations  of  the  form  given  in  (1.7)  for 
the  carrier  transport,  and  the  Helmholtz  eigenvalue  equation  for  the  optical  field.  We 
add  the  heat  flow  equation  given  in  (1.6)  to  obtain  the  temperature  distribution.  We 
henceforth  call  the  first  three  equations  and  the  heat  flow  equation  the  electronic  part 
equations.  Then  the  partial  differential  equations  for  the  electronic  part  are 


-V.eVrJ,  +  q(n-p-N£  +  Ni)  =  0,  (2.1) 

dn/dt  +  V  •  je  +  +  UZ  +  Uhs r  +  t/Aug  =  0,  (2.2) 

dp/dt  +  V  •  jh  +  C&,  +  U&  +  Uhsr  +  ^ug  =  0,  (2.3) 

du/d t  -  V  •  «VT  +  V  •  J E  +  hu  (lEd  -  i?fc)  +  =  0.  (2.4) 


Hore  e  is  the  static  permittivity  of  the  material,  and  N£  |Aa]  is  the  density  of  ionized 
donor  [acceptor]  impurities. 

By  introducing  only  one  temperature  (for  both  the  electron  gas  and  the  crystal  lattice) 
and  only  one  equation  for  the  heat  transfer,  we  have  assumed  that  phonon  scattering 
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leads  quickly  to  a  local  equilibrium.  This  assumption  is  valid  since  our  scale  of  transient 
analysis  is  in  the  order  of  1  ns,  while  the  momentum  relaxation  time  is  in  the  order  of 
1  ps. 

The  Helmholtz  equation  can  be  derived  from  the  Maxwell  equations  by  ignoring 
spatial  derivatives  of  £w,  the  relative  permittivity  at  optical  frequency  u>,  as  follows: 

(V2  +  -02)<f>  =  0,  (2.5) 

where  c  is  the  light  velocity  in  a  vacuum,  and  <f>  and  f3  are  the  optical  field  and  its  prop¬ 
agation  constant,  respectively.  Here  we  solve  the  two-dimensional  Helmholtz  equation 
directly  rather  than  employing  the  effective  refractive  index  method  [28].  By  doing  this 
the  simulator  can  exactly  analyze  virtually  any  refractive  index  profile.  It  should  be 
noted  however  that  this  equation  totally  neglects  the  vectorial  nature  of  the  electromag¬ 
netics  problem.  We  then  simply  use  the  Dirichlet  boundary  condition  at  the  natural 
boundary  assuming  that  the  field  amplitude  at  the  boundary  is  negligible.  When  the 
device  profile  is  symmetric,  we  use  the  vanishing  Neumann  boundary  condition  along  the 
line  of  symmetry,  and  solve  the  equation  in  a  half  domain  for  a  symmetrically  designed 
profile.  We  fix  u  at  the  center  frequency  of  stimulated  emission,  and  solve  (2.5)  for  the 
eigenvalues  (32  and  the  corresponding  eigenfunctions  <f>.  Then  we  get  multiple  transverse 
mode  solutions.  We  cannot  actually  predict  how  those  transverse  modes  will  develop. 
There  are  also  multiple  longitudinal  modes  which  form  a  beat  pattern  in  space,  making 
the  density  of  carriers  irregular  along  the  longitudinal  direction.  In  our  simulation,  we 
have  assumed  that  all  the  characteristics  are  uniform  along  the  optical  axis  either  because 
such  an  irregular  beat  pattern  is  random  so  that  the  overall  effect  can  be  ignored,  or  be¬ 
cause  only  one  optical  mode  is  dominant,  as  in  well-designed  quantum-well  lasers.  Even 
for  a  double-heterojunction  laser,  when  the  laser  is  shorter  than  50  /im  long,  we  need  not 
consider  such  a  problem  since  there  is  only  one  longitudinal  mode  actually  lasing. 
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Outer  iteration  is  done  with  the  photon  rate  equation  in  (1.5),  viz., 


d  S„/dt  =  GUSV  +  Ef  —  Sv(tv. 

We  have  used  a  form  which  explicitly  involves  the  mode  occupation  number  Sv  of  each 
mode  v.  In  principle,  we  can  get  a  self-consistent  solution  if  we  solve  the  above  six 
equations  iteratively. 

2.3  Modeling  of  physical  Parameters 

2.3.1  Sources  of  Poisson  equation 

Throughout  the  formulation,  we  intentionally  use  the  variables  tjc  and  r/h  rather  than 
carrier  densities,  n  and  p,  which  can  be  expressed  in  terms  of  the  former  variables  as  in 
(1.8).  These  expressions  are  valid  in  the  regions  whose  band  edges  can  be  approximated 
to  be  parabolic.  However,  in  size-quantized  regions,  i.e.,  inside  the  quantum  well,  we 
need  to  use  (1.22)-(1.23)  with  m ^  =  0.45r»o  and  m jjf,  =  0.08mo  according  to  [29]  for 
GaAs. 

In  the  semiconductor  laser  application,  we  cannot  ignore  incomplete  ionization  of 
dopants.  Specifically  we  use,  for  N£  and  , 

i\i+  — _ _  ]\i-  = _ _  (9  61 

D  1  +  2exp(?7e  +  eD/T)’  a  1 -f  4exp(%  +  £a /T)’ 

where  cd  and  are  the  absolute  values  of  the  energy  levels  of  donor  and  acceptor  states, 
respectively,  with  respect  to  the  corresponding  band  edges. 
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2.3.2  Hall-Shockley-Reed  recombination 

For  the  modeling  of  the  Hall-Shockley-Reed  (HSR)  recombination  rate,  we  use 

Uhsr  =  — ,  — r,  (2.7) 

n0  =  N'cT^Fin({F-Ec)  /T), 

Po  =  KT3I2?1/2((Ev  -  T)  /T), 
n,  =  NiT^/MEi  -  Ec)/T), 

Pi  =  KT2'2E1/2([EV  -  Et)/T), 

?s(f.+  H)/ 2 

where  n0po  is  the  product  of  equilibrium  densities  of  electrons  and  holes,  while  nx  and  pi 
are  the  respective  densities  when  the  Fermi  energy  is  positioned  at  the  trap  energy  level 
Et.  The  HSR  recombination  is  a  typical  nonradiative  recombination,  and  this  generates 
heat  when  the  net  effect  is  a  recombination  of  electron-hole  pairs.  The  strength  of  this 
is  modeled  by  the  parameters  t„  and  rp  in  (2.7).  The  numerical  values  for  this  depends 
on  the  purity  of  the  device  and  crystal  defects.  In  the  simulation  example,  we  use 
r„  ~  rp  ~  5  •  10-8  s  throughout  the  whole  device,  assuming  that  the  device  is  free  of 
defects  so  that  the  HSR  recombination  is  not  a  major  recombination  process.  Obviously, 
the  values  for  r„  and  rp  can  be  changed  locally  if  necessary. 

2.3.3  Auger  recombination 

It  is  known  that  Auger  recombination  is  not  as  crucial  for  the  temperature  dependence 
of  the  threshold  current  in  GaAs-AlGaAs  lasers  as  it  is  in  long-wavelength  lasers  (e.g., 
InGaAsP  lasers).  There  are  several  different  kinds  of  Auger  processes  according  to  which 
bands  are  involved — CCCH,  CHHS,  and  CHHL  processes — and  to  whether  the  phonon  is 
involved — phonon-less  CHHS  and  phonon-assisted  CHHS  processes.  For  semiconductors 
with  Eq  >  A,  the  CHHL  process  is  negligible  compared  to  the  CHHS  process.  For 
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GaAs,  the  Auger  coefficient  for  phonon-less  CCCH  process  is  zero,  since  the  condition  of 
energy-momentum  conservation  is  not  satisfied  [30],  [31].  Among  the  various  kinds  of  pho¬ 
non-scattering  which  assist  the  Auger  processes,  the  deformation-potential  LO-phonon- 
scattering  is  known  to  be  most  responsible  for  Auger  recombination  at  low  temperatures, 
followed  by  the  polar-optical  LO-phonon-scattering  [32],  [33].  These  two  kinds  of  phonon- 
assisted  Auger  processes  yield  relatively  weak  temperature  dependence.  Consequently, 
we  only  have  the  LO-phonon-assisted  CCCH  process  for  Cn,  whereas  we  have  the  phonon¬ 
less  and  the  LO-phonon-assisted  CHHS  processes  for  C ^  in  the  following  model: 

CJ*U‘=  [C„p+C;hnr](nrp-nj;po),  (2.8) 


where  nr  is  the  electron  density  in  the  T  valley. 

For  ALcGa^As,  we  take 

c,(*)  =  Cp(0)  { [£S(0)  -  A(0)]  /  [£$(*)  -  A(x)]}’A  (2.9) 

C„(x)  =  C„(0)  [£a(0)/£S(x)P  (2.10) 

with  7a  a  1.5  from  approximate  expressions  in  [31]. 

For  the  quantum-well  active  layer,  the  work  of  Takeshima  [34]  suggests  that  C„  <C  Cp 
also  for  the  quasi-two-dimensional  gases.  Since  the  quantum-well  active  layer  is  normally 
undoped,  and  has  almost  the  same  densities  of  electrons  and  holes,  we  thus  can  ignore  the 
CCCH  process.  It  was  found  that  the  phonon-assisted  CHHS  process  for  a  quantum-well 
laser  has  only  a  weak  temperature  dependence  [34].  The  well-thickness  dependence  of 
Auger  processes  at  room  temperatures  was  also  found  to  be  weak  around  the  threshold 
injection  carrier  densities  of  a  quantum-well  laser  [35],  [34].  Hence,  for  the  GaAs  quan¬ 
tum-well  active  layer,  we  use  the  room-temperature  Auger  coefficients  of  the  bulk  GaAs 
for  all  temperatures  for  Cp  in  the  first  approximation. 
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Table  2.1:  GaAs  Characteristics 


parameters 

symb. 

data 

ref. 

Spontaneous  recomb,  coeff. 

Bo 

0.77  •  10-10  cm3/s 

[36] 

CCCH  Auger  coeff.  at  300  K 

Cn 

10'36  cm6-s-1 

[37] 

Phonon-less  CHHS  Auger  coeff.  at  600  K 

c*ph 

2.  •  10~29  cm6  •  s"1 

[38] 

Phonon-assisted  CHHS  Auger  coeff.  at  100  K 

cf 

10-32  cm6  •  s-1 

[39] 

Energy  band  gap  at  0  K 

Eg(  0) 

1.519  eV 

[40] 

Energy  bandgap  temperature  parameters 

aEa 

5.405-4  eV  •  K"1 

[40] 

0 

204  K 

[40] 

Temper,  coeff.  for  free  electron  absorption 

7* 

0.34 

[41] 

LO  phonon  energy 

Threshold  momentum  parameter 

hJ? 

29.6  meV 

[39] 

for  CHHS  Auger  process 

a 

0.22 

[31] 

[100]-plane  quantum-well  eff.  heavy  hole  mass 

mhh 

0.45mo 

[29] 

[100]-plane  quantum- well  eff.  light  hole  mass 
Electron  mobility  parameters 

mlh 

0.08mo 

[29] 

for  ionized-impurity  scattering 

Njon 

9.85  •  1016  cm"3 

[42] 

a*on 

0.553 

[42] 

Hole  mobility  at  300  K  [cm2 /V  s] 

Hole  mobility  param. 

/*h 

400 

[43] 

for  ionized-impurity  scattering 

Nr 

6.25  •  1017  cm"3 

[44] 

Carrier-density  dependence  of  refractive  index 

dy/np 

-1.05  •  10-2°  cm3 

[45] 

Thermal  coeff.  of  refractive  index 

4.9  •  10~4  K'1 

[46] 

Free  carrier  absorption  cross  section 

<tC 

3  ■  10-18  cm2 

[47] 

7  •  10~18  cm2 

[47] 

The  physical  parameters  used  for  the  computation  of  Auger  recombination  for  GaAs 
are  shown  in  Table  2.1. 


2.3.4  Radiative  recombination 

Above  threshold,  we  have  stimulated  emission  as  well  as  spontaneous  emission.  Either 
type  of  light  emission  comes  from  both  the  heavy-hole  ( v  =  e-hh)  and  light-hole  ( v  =  e-lh) 
transitions.  Under  the  assumption  of  strict  fc-selection,  the  mode  gain  and  spontaneous 
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emission  factor  in  (1.5)  can  be  expressed  as 


Gv 


R? 


-i 


(2.11) 


Jv\*,fd3 


{* 

•  =  £  9,(K,h  -  Kj.)  • 

heu  - 

r  i  ( 

KPJ 

v=e-fah,  e-lh 

1 -/(£!»). 

(2.12) 


B„  = 


m\e  quJvM2 


iA/br, 


|  Mb  |  ~ 


m’Eg  1  +  A/^ 


(2.13) 


12m£  1  +  2A/3 

where  i?„  is  the  Einstein  coefficient  of  optical  transition  for  mode  i/,  and  Af  is  the  refractive 
index  around  E  =  %uv.  Integration  denoted  by  fv  is  done  over  the  whole  volume  of  the 
laser  cavity.  The  function  f(E)  is  the  probability  of  an  electron  state  at  energy  E  being 
occupied,  and  the  constant  |Mb|2  is  the  mean  absolute  square  of  the  matrix  element 
of  the  momentum  (to  the  optical  polarization  direction)  obtained  with  respect  to  the 
Bloch  wavefunctions  for  the  conduction  band  electron  state  and  a  randomly  polarized 
valence  band  electron  state  both  at  He  =  0  [48]  [49].  We  normally  use  the  Fermi-Dirac 
probability  function  for  /(•)  by  assuming  that  the  local  equilibrium  is  always  maintained. 
The  function  gv(E*<v  —  Eyv)  (with  v  —  e-hh,  e-lh)  is  the  reduced  density  of  states  for 
optical  transition  between  two  levels  of  energy — E1V  in  the  conduction  band  and  Eyv  v  in 
the  valence  band. 


Kf.^  =  £5 


+ 


~  Eg  nv 

1  +  mr/mhb’  ^'e‘hh 


=  Ey  - 


lxu)i/  —  -F'g 
mhh/ml  +  1 ' 


The  energies  ££e_ ^  and  Ele^  are  similarly  defined  with  mhh  replaced  with  mu,  in  the 
above  equations.  We  optionally  use  a  staircase  quasi-two-dimensional  density  of  states 
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for  the  quantum-well  layer  [50]  and  a  square-root  density  of  states  for  the  other  regions. 


gv(E)dE  =  < 


Ny(E)  d E,  in  quantum-well, 

— — ^ \\j2mv  (E  —  Eg)  dE,  otherwise, 
ir2n 


N;(E)  at  (£-££))  , 

for  t;  =  e-hh,  e-lh  with 


(2.14) 


(2.15) 


”Vhh  =  (l/ml  +  1/mu,)  ‘ ,  me.u,  a  (l/ra[  +  l/mlk)  , 

">’hh  =  (l/m^  +  1/mJ,)  ' ,  m%m  5  (l/m^  +  1/mJ,)  ' , 

where  fq  is  the  thickness  of  the  quantum-well  layer.  The  effective  masses  and 
are  the  heavy  and  light  hole  masses,  respectively,  in  the  quantum-well  layer  parallel  to 
the  [100]  plane.  Function  Int(x)  takes  the  maximum  non-negative  integer  less  than  x. 

Experimental  results  do  not  show  strict  fe-selection.  The  opposite  end  to  this  model 
is  no  Jfe-selection  except  the  subband  selection,  whose  gain  coefficient  spectrum  shape  was 
found  to  agree  better  with  experiment  than  strict  Ks-selection  in  the  case  of  quantum-well 
lasers  [51].  However,  no-fc-selection  model  introduces  several  fitting  parameters  which 
are  multiplication  factors  for  the  momentum  matrix  elements  [52],  [51].  We  use  a  relaxed 
^-selection  model  based  on  recent  publications  [53],  [54],  which  will  need  an  elaborate 
verification  in  the  future.  To  simplify  the  convolution  integral  which  is  not  acceptable 
in  this  two-dimensional  simulation,  we  use  the  following  multiplication  factor  pTclx(E) 
which  is  multiplied  to  the  reduced-density  of  states  in  (2.14).  That  is,  for  the  transition 
between  the  conduction  and  the  heavy  hole  bands, 


1  +  tanh 


( 


(2.16) 
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where  is  a  parametric  relaxation  energy  for  our  approximation.  We  have  chosen 
Enix  —  5  meV  in  order  to  fit  the  results  of  [53].  Note  that  this  relaxation  reduces  the 
maximum  gain  coefficient  substantially,  which  in  turn  increases  the  threshold  current 
estimation.  For  v  =  e-lh  of  light-hole  transition,  £££  ,,  and  are  replaced  by 
E^,  ,  and  m^,  respectively. 

We  have  also  incorporated  polarization-dependent  gain  coefficient  change  for  quan¬ 
tum-well  lasers  according  to  [55].  Between  the  two  orthogonal  polarization  modes  of 
guided  optical  fields,  the  transverse  electric  (TE)  mode  is  usually  dominant  either  because 
of  the  reflectivity  difference  from  cleaved  facets  without  special  coating  [56]  or  because 
of  the  low  chance  of  the  gain  spectrum  over  the  first  light  hole  band  being  the  peak.  We 
therefore  confine  the  treatment  to  that  of  the  TE  mode  transition.  For  the  transition 
between  the  conduction  and  the  heavy  hole  bands, 

,»(*)-!  (l  +  ®±||i).  (2.1T) 

This  extra  factor  is  then  multiplied  to  the  reduced  density  of  states  as  well  as  pnl*(E). 
For  the  transition  between  the  conduction  and  the  light  hole  bands, 

,™(E)  =  !( 2-^l^i).  (2.18) 

The  temperature  dependence  of  the  energy  band  gap  of  the  active  region  material  is 
of  importance  since  the  band  gap  determines  the  wavelength  of  the  stimulated  emission. 
We  take  Varshni’s  [57]  formula: 

Ea(T)  =  Ea(0  K)  -  aEoT2/  (T  +  0) .  (2.19) 


The  radiative  recombination  rate  can  also  be  represented  in  terms  of  r®‘  and  r®p  as 


Ikd  =  E  14  +  /  r r  do,. 


(2.20) 
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The  experimental  data  for  the  spontaneous  radiative  recombination  factor  B^(x)  are 
not  available  for  AlrGaj_r  As.  However,  we  can  deduce  that  it  has  only  a  small  variation 
which  comes  from  the  the  functional  dependence  of  the  Einstein  A  coefficient.  That  is, 

=  BZ(0)  [£a(x)|Mb(x)|2/££(0)|Afb(0)|!]  ,  (2.21) 

where  |JV/b(x)|2  is  given  in  (2.13)  with  ar-dependent  and  Eq. 

The  numerical  integration  in  (2.20)  imposes  a  heavy  computational  load.  We  therefore 
compute  the  second  term  of  the  right-hand  side  of  (2.20)  with  (2.12)  inside  the  active 
layer  only,  where  most  radiative  recombination  takes  place.  In  the  other  regions,  we 
approximate  (2.20)  into 

Urad  =  +  B0  ( nrp  -  n£po)  • 

U 

The  value  of  the  constant  B0  is  given  in  Table  2.1. 


2.3.5  Refractive  index 

To  accommodate  the  refractive  index  variations  for  all  the  regions  including  the  active 
layer,  we  use  the  following  formula: 


(F  t\\T  xk-1"*5  2 

t-=)  v/Sp  +  ^T  (T-7i)  +A4J, 

“v”P/r=T0 


(2.22) 


a  -fc  92  ( nT  n X  nL  P  \ 

— - J  I  f  +  — T  - L  - I  ’ 

EqU)2  \m‘  m*  m^j 


where  e0  is  the  vacuum  permittivity,  {3?  is  the  thermal  coefficient  of  refractive  index 
change,  and  (dA/’ /dy/™p)T_T  represents  the  change  due  to  band-to-band  transitions  of 
carriers  via  the  Kramers-Kronig  relation.  The  latter  value  for  the  active  region  is  obtained 
by  subtracting  the  theoretical  free-carrier  components  from  the  empirical  data  of  [45]  for 
the  overall  carrier-density  dependence.  In  (2.22),  v/np  has  been  taken  to  interpret  the 
empirical  values  obtained  under  the  approximation  n  ~  p  in  the  active  layer.  For  the 
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other  regions,  we  neglect  the  change  in  (2.22)  due  to  the  band-to-band  transition,  either 
because  the  absorption  edge  is  well  above  the  stimulated  emission  frequency  or  because 
the  optical  field  intensity  is  negligibly  small  in  those  regions. 

For  strongly  gain-guided  structures,  we  should  include  change  of  the  imaginary  part 
of  the  refractive  index.  This  results  in  a  complex  matrix  of  a  non-Hermitian  type  after 
discretizing  the  Helmholtz  equation,  (2.5).  The  work  for  the  gain-guided  structures  is 
reserved  for  future  research. 


2.3.6  Free-Carrier  absorption  and  photon  lifetime 

The  free-carrier  absorption  rate  is  proportional  to  the  light  intensity  and  the  carrier 
densities.  It  is  related  to  the  free-carrier- absorption  coefficient  afc,  evaluated  for  a  balk 
semiconductor,  as  follows: 

Rtc  =  -^raic  Sv  \j>*\2 ,  a{c  =  akn  +  akp, 

•'’’g  i/ 

where  a and  er£c  is  the  free-carrier-absorption  cross-sections  for  electrons  and  holes, 
respectively.  The  temperature  variation  of  these  coefficients  for  bulk  GaAs  is  not  strong. 
For  instance,  in  the  case  of  electrons,  we  can  only  deduce  r£c  ~  T0  3  around  300  K  from  the 
data  in  [41].  We  therefore  ignore  temperature  dependence  of  the  free-carrier-absorption 
coefficients,  and  take  a  simple  model  [47]  of  constant  <r„  and  ak. 

For  the  graded-j unction  cladding  layers  which  sandwich  the  active  layer,  to  our  knowl¬ 
edge,  experimental  data  are  not  yet  available  for  the  various  mole  fraction  x  of  AlAs.  We 
thus  proceed  with  approximate  values  deduced  from  the  available  GaAs  data. 

The  overall  free-carrier-absorption  coefficient  for  the  calculation  of  photon  lifetime  rv 
follows  from  the  above  absorption  rate  as 

a'c  =  fvakm2d3r.  (2.23) 
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The  overall  inverse  photon  lifetime  which  we  need  for  (1.5)  is  then  expressed  as 

TV-1  =  v*  Kc  +  <C8  +  (2L)-1  ln(l /R,R2)\  ,  (2.24) 

where  v6  is  the  group  velocity  of  mode  i/,  and  a*‘g  is  the  waveguide-scattering  loss 
coefficient  for  mode  u,  which  depends  on  geometrical  irregularity  formed  during  the 
fabrication  process  of  the  laser  waveguide. 

Since  there  are  graded  heterojunction  regions  in  the  profile  of  semiconductor  lasers,  we 
need  values  of  material  parameters  as  functions  of  variable  A1  composition.  In  Table  2.2, 
such  parameters  which  are  obtained  with  a  linear  or  quadratic  interpolation  are  shown.  In 
the  table  the  lattice  thermal  conductivity  k  in  (2.4)  is  represented  by  the  inverse,  whose 
quadratic  approximation  fits  better  with  experiments  than  that  of  k  itself.  However, 
parameters  such  as  carrier  mobilities  are  not  as  simple  to  be  represented  as  others  in 
Table  2.2. 


2.3.7  Carrier  mobilities  of  AlIGa1_rAs 

The  electron  mobility  reflects  the  lowest  band  minimum  among  three  conduction  band 
minima,  T,  L,  and  X.  We  seek  to  express  the  experimental  results  in  [68]  as  follows: 


8000  -  1.818  •  104x,  0  <  x  <  0.429, 

90  +  1.145  •  10s  (x  -  0.46)2 ,  0.429  <  x  <  0.46, 

90  +  3.75  •  104  (x  -  0.46)2 ,  0.46  <  x  <  0.50, 

200  -  2.0/ (x  -  0.46) ,  0.50  <  x  <  1, 


(2.25) 


in  units  of  cm2/V-s.  The  temperature  dependence  of  electron  mobility  in  doped  AlGaAs 
(or  even  GaAs)  is  not  well  modeled  due  to  various  factors  including  the  above-mentioned 
complex  conduction  band  structure  [63].  The  modeling  of  this  is  reserved  for  future  work. 

scattering-limited  Here  we  multiply  a  crude  doping-dependence  factor  to  the  expres¬ 
sion  for  the  electron  mobility  of  lightly  doped  AlGaAs  in  (2.25).  This  factor  was  taken 
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Table  2.2:  AlxGai_rAs  Characteristics 


parameters 

symb. 

function  of  z,  A1  composition 

ref. 

Static  permittivity 

i 

13.18 +  (10.06-  13.18)  x 

[58], [59] 

Refractive  index  at  ku 7^  =  1.38  eV 

AT 

3.590  -  0.710x  +  0.091x2 

[47] 

High  freq.  permittivity 

€qo 

10.89  +  (8.16  -  10.89)  x 

[58], [59] 

T  valley  enegy  gap  at  300  K  [eV] 

Erc 

1.424  +  (2.671  -  1.424)  x 

(47! 

and  its  temper,  change  [eV/K] 

d Sfi 
dT 

[-3.95  -(5.1  -  3.95)  x]  •  10"4 

[60], [61] 

X  valley  energy  gap  at  300  K  [eV] 

Eg 

1.9 +  (2.168  -  1.9)  x 

[62], [63] 

and  its  temper,  change  [eV  /K] 

Diff.  between  the  valence  band 

d E% 
d  T 

-3.6  •  10-4 

[61], [63] 

and  spin-split  band  edges  [eV] 
Donor  ionization  energy 

A 

.340  -  .040x 

[63] 

for  T  band  [meV] 

Donor  ionization  energy 

er 

5  +  5x 

[64], [63] 

for  X  band  [meV] 

e* 

25  +  lOx 

[64], [63] 

LO  phonon  energy  [meV] 

36.25  -  6.55x  +  1.79x2 

[65] 

T  valley  electron  mass 

X  valley  density-of-states-equiv. 

< 

.067 +  (.15  -. 067) x  [m0] 

[47] 

electron  mass 

m* 

.32  +  (.26  —  .32)  x  [m0] 

[63] 

Heavy  hole  mass 

.62  +  (.76  -  .62)  x  [m0] 

[66] 

Light  hole  mass 

™lh 

.087 +  (.15  -. 087) x  [m0] 

[66], [67] 

Thermal  resistivity  [K-cm/W] 

w 

2.27  +  28.83x  -  30x2 

[63] 
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from  data  in  [42]  for  GaAs., 


M,  = 


0.34  + 


0.66 


1  + 


( 


N$  +  Nj  +PhhV!!< 

Nr  ) 


(2.26) 


where  phh  is  the  density  of  heavy  holes  which  act  like  ionized  impurities  for  electrons  [69]. 

To  our  knowledge,  no  experimental  data  published  for  the  compositional  variation  of 
the  hole  mobility  in  AlrGai_xAs  are  available.  We  generalize  the  theoretical  scheme  in 
[70]  for  the  mobility  parameter  Mh  as 


1 

1 

-l/£s(x)  -  l/e0O(x)  , 

(Li 

Mh"' 

x  =  0 

[l/e8(0)  —  l/£oo(0)  1 

l To ) 

T  =  T0 


jVp  +  Na  mhh(x)  +  mih(g)  ( g,t(0)  V  (To'? 
Ahon(0)  mhh(0)  +  mih(O)  V£st(x)/  \T  ) 


(2.27) 


where  N^a  refers  to  the  parameter  for  the  ionized-impurity-scattering-limited  hole  mo¬ 
bility  for  GaAs,  and  mhh(x)  and  mu,( x)  represent  the  appropriate  effective  masses  with 
the  A1  composition  x.  £<»  is  the  dielectric  constant  at  high  frequency,  whose  normalized 
value  is  found  in  Table  2.2. 


2.4  Numerical  Approach 

The  box-discretization  [20]  of  these  device  equations  has  been  used  with  the  gen¬ 
eralized  Scharfetter-Gummel  discretization  scheme  [21]  described  in  Sec.  1.5.  The  four 
equations  (2. 1 )— (2.4)  for  the  electronic  part  are  then  solved  by  the  full-Newton  method 
and  a  sparse  matrix  package  [71]— [72].  Since  the  current  flow  near  the  surface  does  not 
count  for  the  overall  characteristics  in  the  case  of  semiconductor  lasers,  simple  ohmic 
contacts  at  the  electrodes  have  been  assumed,  and  Neumann  boundary  conditions  with 
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vanishing  flux  have  been  used  at  the  surfaces  without  electrode  contact  throughout  the 
simulation.  Current  boundary  conditions  and  the  heat  current  boundary  conditions  have 
been  employed  in  consideration  of  application  practice  in  driver  circuits  for  semiconductor 
lasers. 

The  steady-state  solution  can  be  obtained  by  solving  the  steady-state  problem  for  the 
electronic  part  equations  by  letting  d/d/  zero  and  solving  the  time-dependent  problem 
for  the  rate  equation.  Other  authors  of  semiconductor  laser  simulations  used  a  simpler 
method  which  was  described  in  [16].  However,  when  there  axe  multiple  modes  competing 
with  each  other,  such  methods  cannot  be  used. 

Since  we  have  adopted  the  full-Newton  method  to  solve  the  electronic  part  equations, 
we  can  employ  the  fully-implicit-backward  Euler  method  [14]  to  obtain  the  dynamic 
solver,  thus  avoiding  the  stability  problem  associated  with  the  first-order  time  derivatives. 
The  rate  equation  is  also  treated  with  the  full-backward  Euler  method,  and  this  rate 
equation  is  coupled  with  the  electronic  part  equation  by  solving  the  two  time  dependent 
problems  alternatively.  The  optical  field  pattern  is  obtained  by  solving  the  Helmholtz 
eigenvalue  equation  at  each  iteration.  We  should  understand  that  the  whole  procedure  of 
time  dependent  solution  then  does  not  become  a  fully  implicit  backward- Euler  procedure 
since  we  still  have  two  parts  to  be  solved  iteratively.  This  limits  the  size  of  the  time  step 
which  is  set  especially  during  the  transients  of  the  relaxation  oscillation  which  appears 
when  the  laser  is  driven  by  a  step  current. 

The  Helmholtz  eigenvalue  equation,  (2.5),  is  solved  by  Rutishauser’s  subspace  itera¬ 
tion  subroutine  RITZIT  [73].  The  routine  solves  for  the  desired  number  of  eigenvalues  of 
the  largest  magnitude  and  the  corresponding  eigenvectors.  From  the  mathematical  na¬ 
ture  of  the  problem,  the  largest  magnitude  eigenvalues  are  negative.  Therefore,  we  need 
to  shift  the  spectrum  to  make  the  most  positive  eigenvalue — the  fundamental  transverse 
mode  eigenvalue — the  dominant  one  [74].  We  then  need  not  employ  the  inverse  itera¬ 
tion.  This  method  allows  us  to  obtain  only  the  desired  physical  eigenvalues  and  their 
characteristic  optical  field  patterns. 
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The  overall  flowchart  of  the  program  is  shown  in  Figure  2.1.  It  has  an  input  card  pars¬ 
ing  structure.  Graphical  outputs  are  generated  exclusively  in  the  PostScript  language 
source  code  format,  which  allows  a  flexible,  device-independent,  and  unilateral  operation 
for  all  types  of  graphical  output. 


2.5  Application  to  Quantum- Well  Lasers 

We  have  applied  our  simulator  to  a  380  /xm-long  GaAs-AlGaAs  GRINSCH  buried- 
quantum-well  laser  [27],  wnose  cross-sectional  view  is  shown  in  Figure  2.2.  The  semi¬ 
conductor  equations  yield  the  three  solutions  for  the  potential,  the  electron  and  the  hole 
densities.  The  temperature  distribution  over  the  cross-sectional  profile  is  depicted  in 
Figure  2.3.  We  have  assumed  that  the  cathode  is  in  contact  with  an  ideal  heat  sink,  and 
that  the  energy  flow  out  of  the  anode  side  is  negligible.  We  have  found  that  electrons 
and  holes  are  very  effective  energy  carriers.  The  increase  of  the  temperature  slope  in 
the  middle  of  the  cross  section  represents  the  poor  thermal  conductivity  of  AlGaAs  com¬ 
pared  to  pure  GaAs  in  the  substrate.  The  hump  above  the  quantum-well  active  region 
represents  a  heat  dissipation  from  the  electron  system  by  the  free-carrier  absorption  of 
the  photon  energy.  The  temperature  elevation  around  the  active  layer  depends  strongly 
on  the  distance  between  the  active  layer  and  the  heat  sink. 

The  current  flow  diagram  is  shown  in  Figure  2.4.  The  figure  shows  that  the  re¬ 
verse-biased  junction  across  the  blocking  region  of  higher  bandgap  sufficiently  channels 
the  current  through  the  narrow  forward- biased  region  without  introducing  any  insula¬ 
tor  material  into  the  laser  profile.  The  thermionic  injection  current  simulation  has  been 
found  indispensable  in  the  simulation  of  a  quantum-well  laser,  since  the  different  treat¬ 
ment  of  normalized  carriers  in  the  discretized  carrier  transport  equations  as  described  in 
(1.30)  in  the  different  regions  of  the  device — inside  and  outside  the  quantum  well — yields 
an  unphysical  result  with  only  a  drift-diffusion  model  applied  everywhere.  Figure  2.5 
demonstrates  this  by  showing  the  two  uninterpolated  pictures  of  current  flow  diagrams, 
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Figure  2.1:  Flowchart  of  MINILASE. 


Figure  2.2:  (a)  Schematic  cross-section  of  a  simplified  GaAs-AlGaAs  graded-index- 
separate-confinement-heterojunction  buried-quantum-well  laser.  The  layers  above  and 
below  the  active  layer  are  quadratically  graded  from  Al.25Ga.75As  to  Al.5Ga.5As.  These 
three  layers  are  undoped,  (b)  Variation  of  A1  composition  along  the  center  line. 

one  with  a  thermionic  emission  simulation  and  the  other  with  a  drift-diffusion  model 
only. 

We  have  also  obtained  a  light-current  curve  as  shown  in  Figure  2.6  by  varying  the 
carrier  injection.  This  shows  a  threshold  current  of  ~  3.75  mA  and  a  local  threshold 
current  density  of  380  A/cm2  inside  the  channel  between  the  two  blocking  regions  of  the 
higher  bandgap.  The  latter  number  matches  very  well  with  various  experimental  data 
of  broad-area  single-quantum-well  GRINSCH  lasers  [75]— [76].  The  values  from  others’ 
theoretical  calculations  [77],  [50]  vary  from  240  A/cm2  to  ~  550  A/cm2,  since  each  author 
used  different  values  for  the  material  constants  (e.g.,  the  waveguide  scattering  loss 
coefficient,  and  the  intrinsic  loss  coefficient,  and  also  treated  the  Jfe-selection  rule  and 
polarization  dependence  differently. 

From  the  eigenvalue  solver,  we  have  obtained  the  optical  field  profile  as  shown  in 
Figure  2.7.  From  the  spectral  analysis,  we  have  obtained  the  gain  spectra  (Figure  2.8) 
and  the  optical  output  spectra  (Figure  2.10)  for  various  carrier  injection  levels.  The 
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temperature  [K] 


Figure  2.3:  Temperature  distribution  profile  for  the  half  domain  at  a  180  mA  current 
level.  The  cathode  is  assumed  to  be  in  contact  with  an  ideal  heat  sink.  Energy  flow  out 
of  the  anode  side  is  assumed  to  be  negligible. 
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Figure  2.4:  Current  flow  diagram  for  the  half  domain  at  4.0  mA.  One  fortieth  of  the 
total  current  (electron  current  plus  hole  current)  flows  between  the  two  adjacent  lines  in 
the  diagram. 
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Figure  2.5:  Uninterpolated  current  flow  diagram  for  the  half  domain  at  3  mA  with 
simulations  (a)  of  thermionic  emission  as  well  as  drift-diffusion  and  (b)  of  drift-diffusion 
only. 
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Figure  2.6:  Light-current  curve  of  a  model  quantum-well  laser,  (a)  Around  the  thresh¬ 
old.  (b)  Over  the  injection  current  of  0-160  mA. 
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Figure  2.7:  Optical  intensity  profile  around  the  quantum-well  region  yielding  an  optical 
confinement  factor  of  3.9  %. 
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Figure  2.8:  Gain  coefficients  at  various  current  levels.  Free-carrier  loss  has  been  ob¬ 
tained  to  be  2.0  /cm,  and  the  waveguide  scattering  loss  has  been  set  at  5.0  /cm.  (a) 
Around  threshold  level.  The  results  show  the  gain  saturation  above  the  threshold. 

(continued  on  p.  41) 
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QW  thickness  20  nm,  laser  length  0.38  mm 


Figure  2.8:  (continued)  (b)  Far  above  threshold.  The  gain  curve  shifts  to  the  lower 
energy  spectrum  as  the  current  increases,  because  the  band  gap  changes  due  to  increasing 
temperature. 


well  temperature  [K] 


gain  coefficient  (G„/ug  with  G„  in  (2.11)  and  the  group  velocity  vs)  curves  show  the  gain 
saturation  effect  above  the  threshold.  Also,  the  gain  curve  shifts  to  the  lower  energy 
spectrum  as  the  current  increases  because  the  band  gap  changes  due  to  increasing  tem¬ 
perature  at  the  quantum  well.  The  temperature  in  the  middle  of  the  quantum  well  is 
shown  in  Figure  2.9.  The  free  carrier  loss  is  obtained  from  a  numerical  integration  as  in 
(2.23),  and  the  value  is  about  2.0  cm-1  above  threshold.  The  waveguide  scattering  loss 
is  still  an  unknown  fitting  parameter.  We  have  chosen  a  typical  value,  5.0  cm"1,  from 
literature  [47]. 

Users  of  MINILASE  can  optionally  choose  the  staircase-like  density  of  states  for  the 
region  assigned  to  be  a  quantum  well  with  quasi-two-dimensional  electron  and  hole  gases. 
The  resulting  gain  and  mode  spectra  from  the  staircase-like  density  of  states  are  shown 
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light  intensity  (occup.  No.) 


Figure  2.10:  Spectra  of  photon  numbers  of  fundamental-longitudinal  modes  at  various 
pumping  levels,  A  staircase-like  density  of  states  for  a  quasi-two-dimensional  electron  gas 
is  used,  (a)  2  mA. 

(continued  on  p.  44) 

in  Figures  2.8  and  2.10,  respectively.  These  diagrams  show  size-quantization  effects  and 
the  subband  level  which  supports  the  stimulated  emission. 

A  mode  development  diagram  can  be  obtained  from  the  transient  simulation,  and  is 
shown  in  Figure  2.11.  We  have  assumed  that  the  laser  is  driven  from  2  mA  abruptly 
to  10  mA  at  t  =  0  ns  and  from  10  mA  back  to  2  mA  at  t  —  7.5  ns.  The  size  of  the 
time  step  has  been  set  uniformly  at  25  ps,  which  has  been  found  to  be  adequately  short. 
The  figure  shows  the  initial  turn-on  delay  time  of  0.9  ns  and  relaxation  oscillations  of 
the  optical  power  for  2  ns.  Also  it  shows  that  several  modes  develop  their  magnitudes 
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Figure  2.11:  Both  turn-on  and  turn-off  transients  in  time.  A  series  of  current  level 
change  from  2  mA  to  10  mA  at  t  =  0  ns  and  from  10  raA  back  to  2  mA  at  t  =  7.5  ns  has 
been  simulated.  Each  line  represents  intensity  development  of  each  mode.  The  central 
mode  is  denoted  by  half  tone  over  the  entire  time  scale  to  show  its  location  in  spectrum 
at  the  initial  stage.  The  time  step  size  of  transient  simulation  has  been  set  at  25  ps. 
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simultaneously  at  first,  but  only  one  mode  has  significant  power  in  steady  state.  Also 
the  center  frequency  and  the  spectral  width  change  during  the  time  development.  The 
turn-on  delay  time  and  the  period  of  relaxation  oscillation  depend  on  the  various  driving 
conditions  of  the  laser.  As  is  evident  from  the  figure,  the  turn-off  transient  is  negligibly 
short  compared  to  the  initial  turn-on  delay  and  the  relaxation  oscillation. 

Figure  2.12  was  obtained  by  applying  a  current  step  from  2  mA  to  50  mA.  The  turn¬ 
on  delay  time  has  been  shortened  to  0.22  ns  and  the  period  of  relaxation  oscillation  has 
become  shorter  than  the  case  of  Figure  2.11.  In  Figure  2.13,  we  show  the  time-develop¬ 
ment  of  carrier  densities  inside  the  quantum  well.  The  magnitude  of  the  initial  peak 
in  the  optical  power  depends  on  the  size  of  the  time  step  taken  for  the  time-dependent 
solution.  This  is  because  the  whole  procedure  of  the  time-dependent  solution  is  still  not 
fully  backward  Euler,  which  was  discussed  in  Sec.  2.4.  Inadequately  long  time  step  gives 
an  excess  peak  in  the  initial  relaxation  oscillation  of  the  optical  power.  For  instance, 
when  we  used  the  time  step  size  of  5  ps  instead  of  2.5  ps  for  the  current  pulse  of  2  -  50  - 
2  mA  sequence,  we  have  observed  the  initial  relaxation  peak  about  10  %  larger  than  the 
peak  observed  with  the  time  step  size  of  2.5  ps.  One  additional  cause  of  the  excessive 
peak  in  the  initial  relaxation  oscillation  could  be  because  we  have  used  a  rate  equation 
which  averages  out  the  intensity  variation  along  the  optical  axis,  and  therefore  ignores 
time  required  to  sweep  through  the  optical  axis. 

All  these  transient  response  characteristics  are  in  good  agreement  with  what  experi¬ 
ments  have  suggested  [78]. 
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Figure  2.12:  Both  turn-on  and  turn-off  transients  in  time.  A  series  of  current  level 
change  from  2  mA  to  50  mA  at  t  =  0  ns  and  from  50  mA  back  to  2  mA  at  t  =  1.5  ns 
has  been  simulated.  Each  line  represents  intensity  development  of  each  mode.  The  mode 
distiguished  from  others  by  half-tone  is  the  one  which  gives  the  peak  at  initial  relaxation 
oscillation.  The  time  step  size  is  set  at  2.5  ps. 
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Figure  2.13:  Development  of  carrier  densities  (n  ~  p)  in  time  in  the  middle  of  the 
quantum  well  in  response  to  square  pulses.  The  change  of  the  laser  light  output  is  also 
shown  with  an  arbitrary  scale,  (a)  In  response  to  a  square  pulse,  2  mA  -  10  mA  -  2  mA. 

(continued  on  p.  50) 
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Figure  2.13:  (continued)  (b)  In  response  to  a  square  pulse,  2  mA  -  50  mA  -  2  mA. 
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CHAPTER  3 


CONCLUSIONS 

We  have  developed  a  two-dimensional  simulator  (MINILASE)  which  analyzes  vir¬ 
tually  all  aspects  of  semiconductor  laser  operations.  The  simulator  solves  the  device 
equations,  heat  flow  equation,  optical  field  equation,  and  the  photon  rate  equation  in  a 
self-consistent  manner.  It  gives  all  the  solutions  to  the  standard  device  equations  and  the 
resulting  current  flow  diagrams  and  the  temperature  profile.  It  also  generates  the  light- 
current  curves,  the  optical  field  profile,  the  temperature  profile,  the  spontaneous-emission 
spectrum,  the  optical  output  spectrum,  and  the  gain  spectrum. 

In  order  to  solve  for  the  internal  temperature  distribution  over  the  device  cross  sec¬ 
tion,  the  energy  transfer  in  degenerate  semiconductor  devices  has  been  examined,  and 
the  relevant  differential  equation  has  been  solved  together  with  the  device  equations. 
Proper  inclusion  of  the  Fermi-Dirac  statistics  has  been  found  to  be  essential  to  obtain 
a  differential  equation  which  is  linearly  independent  of  the  rest  of  the  device  equations. 
A  set  of  novel  expressions  for  the  various  fluxes  has  been  suggested  to  account  for  the 
Fermi-Dirac  statistics  for  the  most  important  scattering  mechanism,  in  which  the  scatter¬ 
ing  time  is  inversely  proportional  to  the  group  velocity  of  carriers.  Explicit  formulas  for 
the  discretization  of  these  fluxes  have  been  found  by  extending  the  Scharfetter-Gummel 
scheme.  There  are  many  limitations  in  practical  circumstances  with  regard  to  the  bound¬ 
ary  condition  of  the  energy  transport  equation.  A  thorough  treatment  for  the  internal 
distribution  of  the  temperature  is  a  complex  problem  in  which  one  needs  to  simulate  the 
outside  of  the  device  as  well  as  the  inside. 

Although  derived  under  a  series  of  assumptions  and  approximations,  we  think  that 
the  formalism  presented  in  this  thesis  is  general  enough  to  be  used  in  simulation  of 
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heat  transfer  in  many  types  of  semiconductor  devices,  including  heterojunction  devices 
such  as  semiconductor  lasers  where  local  power  consumption  is  an  important  factor  for 
device  performance.  The  algorithm  presented  here  has  been  subjected  to  a  series  of 
experiments  on  a  computer  in  the  MINILASE  project  of  the  author,  and  results  were 
given  in  Chapter  2. 

We  have  presented  the  first  two-dimensional  simulation  of  the  optical  processes  in  a 
semiconductor  laser  using  the  rigorous  Einstein  coefficients  of  optical  transition  for  indi¬ 
vidual  Fabry-Perot  modes  of  the  laser  resonator.  Consequently,  it  is  capable  of  simulating 
quantum-well  lasers  and  their  spectral  responses.  MINILASE  is  also  the  first  attempt 
to  use  the  full-Newton  method  for  the  system  of  equations  obtained  from  discretization 
of  the  standard  device  equations  and  the  heat  flow  equation.  This  enables  us  to  analyze 
the  temperature  effect  of  a  forward- biased  semiconductor  device  self-consistently  on  var¬ 
ious  physical  parameters  such  as  the  threshold  current.  Finally,  MINILASE  can  obtain 
the  two-dimensional  transient-response  including  the  initial  relaxation  oscillation  of  the 
individual  Fabry-Perot  mode  intensities  in  a  semiconductor  laser. 

The  first  priority  of  the  MINILASE  program  has  been  a  faithful  implementation  of 
physics,  yet  retaining  versatility  and  flexibility  of  the  program  for  a  wide  range  of  appli¬ 
cations.  We  expect  the  simulator  to  be  a  useful  design  tool  for  optimizing  semiconductor 
lasers  as  well  as  a  useful  research  tool  in  quantum  electronics  laboratories. 
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APPENDIX  A 


ENERGY  CONSERVATION  FROM  THE  BOLTZMANN 

TRANSPORT  EQUATION 


Azoff  [7]  has  derived  the  energy  transport  equation  for  the  nonparabolic  band  struc¬ 
ture  from  the  classical  Boltzmann  transport  equation  (BTE)  by  taking  the  second  mo¬ 
ment  (cited  as  the  third  moment  in  [7]).  The  same  equation  can  be  derived  in  a  much 
simpler  way  by  using  an  equivalent  form  of  the  BTE  in  which  the  velocity  t>k  and  the 
force  Fk  follow  the  divergence  symbols.  Also,  rather  than  taking  the  second  moment 
for  the  standard  equation,  we  can  set  up  a  version  of  BTE  with  respect  to  {E±  -f  Ec)  /k 
for  the  electrons  in  the  conduction  band,  since  the  total  energy  is  conserved  in  much  the 
same  way  as  the  tot’d  number  of  particles  is  conserved.  In  this  way,  no  complication  due 
to  the  nonparabolicity  of  the  band  appears  in  the  derivation  of  the  relation  for  the  first 
law  of  thermodynamics.  That  is, 

^(Ec  +  Fk)/k  +  Vk-^(Fc  +  Fk)/k 

+  V  •  ®k  {E\  -f  Ec)  h  =  7c-  (£k  +  -E'c)/k  •  (A.l) 

idt  c,r 

The  subscripts  c  and  r  stand  for  collision  and  recombination,  respectively.  Integrating 
this  equation  over  the  Brillouin  zone  in  the  reduced-zone  scheme  yields 


I  (£°" + + v  •  (s“" + **) = L E*  (f  )„  S 


-  Ec  (UhSR  +  {/Aug  +  £4a d)  , 


(A-2) 


which  then  reduces  to  (1.1).  Here  tt|fn  and  can  be  identified  as  the  kinetic  energy 
density  and  the  kinetic  energy  flux  of  electrons,  respectively.  Substituting  (1.15)  into 


(A. 2)  yields  the  standard  equation. 


du*n 

dt 


+  V  •  S“n 


+  VEc-je  = 


(A. 3) 


Note  that  the  conserved  quantity  is  not  just  the  internal  energy,  but  the  kinetic  energy 
which  is  the  internal  energy  plus  fsx(Ek  —  Ek_^)f^d3k s/47r3.  Note  also  that  this  expression 
explicitly  separates  out  the  Joule  heat  from  the  total  energy  flux  in  (1.1).  A  source  of 
ill-conditioning  in  possible  implementations  of  this  term  is  in  the  form  of  an  inner  product 
of  the  two  gradients  [79]. 
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APPENDIX  B 


EXPRESSIONS  FOR  THE  HEAT  FLUX 


B.l  From  Relaxation-Time  Approximation 


The  integral  expression  for  the  energy  flux  from  electrons  is 


S'  =  J  vk  (£„  +  Ec)  =  +  Eci.- 


(B.l) 


From  the  relaxation-time  approximation,  the  probability  distribution  function  /k  satisfies 

[12] 

/k  =  (A  +  /-k)  /2,  /k  =  (A  -  /-k)  /2,  (B-2) 


/k  =  ( d/8t  +  «k  •  V  -  VE/h  ■  Vk)  . 


(B.3) 


We  obtain  /k  by  approximately  equating  in  the  right-hand  side  of  (B.3)  to  the  Fermi 
distribution, 

h m  =  {exp[(£k  +  Ec  -  F')  IT']  +  I}"1 , 

where  Tc  is  the  adjusted  electron  temperature  for  nonequilibrium  under  the  influence  of 
an  electric  field  [12].  With  this  /k  ,  we  obtain  the  following  expression  for  Se  in  (B.l)  [5] 


sy-ri,  j.=io, 


I,  =  ~l^kEirk  •  (VF.+  ^VTC) 


(B.4) 

(B.5) 


The  relaxation  times  for  various  scattering  mechanisms  as  functions  of  velocity  or  kinetic 
energy  of  electrons  can  be  found,  e.g.,  in  [8],  [12],  and  [80]. 
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Consider  now  a  case  where  the  relaxation  time  is  proportional  to  jt?k|  1 . 

‘-W-  (B'6) 

where  r  is  velocity-independent.  As  considered  by  Stratton  [1]  with  ®k  =  hk/me,  the 
above  form  is  particularly  suited  for  such  materials  as  moderately-doped  Si  and  Ge,  where 
acoustic  phonon  scattering  represents  a  major  scattering  mechanism.  For  moderately- 
doped  III-V  semiconductors  where  optical  phonon  scattering  is  predominant,  the  relax¬ 
ation-time  approximation  with  (B.6)  is  not  strictly  valid  [12],  [80],  although  still  a  good 
approximation  at  high  temperatures  [12].  Notice  that  ionized  impurity  scattering  fol¬ 
lows  the  same  velocity  dependence  for  the  case  of  heavy  screening,  since  in  this  limit 
the  scattering  potential  is  of  short  range  [12].  Therefore,  for  laser  application,  factoring 
out  l^kl-1  appears  to  be  an  excellent  approximation,  at  least  as  long  as  the  experimen¬ 
tal  mobility  exhibits  degrading  with  increasing  temperature  and  with  increasing  carrier 
densities. 

Following  Kane  [48],  Nag  and  Chakravarti  [81]  derived  an  approximate  relation  be¬ 
tween  k 2  and  Ek  of  an  isotropic  band  structure  under  the  assumption  that  «  £g: 

h2k2/2me  =  Ek(l  +  aeEk ) .  (B.7) 

The  analytic  expression  for  the  nonparabolicity  ae  is  given  in  the  literature  for  various 
semiconductors  [82],  [44],  [83].  Then,  (B.5)  becomes 

=  M.  {[7?  u  +  +  o.r>«  (j  +  2)!^+1(,.)j 

x  (V£c  +  TCV,.)  +  [T>  (j  +  2)!*>+i0fc) 

+  a,T>»  ( j  +  3)LF, „(,«)]  VTt}  .  (B.8) 
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and  the  electron  particle  flux  and  the  energy  flux  are  given  by 


L  =  —Me  {[MV*)  +  2aeJ ^ (i/.)]  (V£c  +  TcVr,e) 

+  [2  Jifoe)  +  6a.*ifo.)]  VTC},  (B.9) 

5e  =  i?cie  -  MeTc  {[^fo.)  +  (V£c  +  TeVi/.) 

+  [6-^2(»?e)  +  24Qe^r3(j;e)]  VTC}  ,  (B.10) 


respectively. 

With  this  expression  for  the  energy  flux,  we  can  separate  the  energy  flux  into  a  com¬ 
ponent  which  is  carried  by  particle  flux  and  the  thermal  conduction  which  is  independent 
of  particle  transfer,  by  rewriting  (B.10)  as 

Se  =  Ecie  +  2r eTcje  ~  ^VTC,  (B.ll) 


_  Ti  +  3Qe^2 
'  “  JF0  +  2 CteK  ' 

Ke  =  MeTcUf2  +  24aejF3  -  -  y  •  (B.12) 

t  f0  +  2are^i  J 

The  thermal  conductivity  Ke  reduces  to  2TC  times  the  electrical  conductivity  for  nonde¬ 
generate  semiconductors  as  expected.  For  extreme  degeneracy  [10], 

fj(v)  —  V3+1/(j  +  1)!  +  7 r2>/-7_1/6j!,  as  r)  — +  oo,  (B.13) 

Ke  reduces  to  7r2Tc/3  times  the  electrical  conductivity,  as  in  the  case  of  a  metal.  Over¬ 
all  the  thermal  conductivity  is  proportional  to  Tc  times  the  electrical  conductivity  in 
agreement  with  the  Wiedemann-Franz  law. 
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Equation  (B.ll)  can  be  used  for  the  simulation  of  semiconductors  neglecting  the  last 
term  when  Ke  is  negligibly  small  compared  to  k  of  the  crystal  lattice.  There  exists 
a  possibility  of  an  ill-conditioned  matrix  problem  when  T  is  very  close  to  1  (as  in  a 
nondegenerate  region  of  a  device)  or  almost  constant  (where  the  grid  structure  is  too 
fine).  Our  conditioning  treatment  in  Sec.  1.4  should  be  helpful  in  those  cases. 

As  was  pointed  out  in  Ref.  [13,  p.  26],  the  term  “heat  flux”  can  be  defined  in  various 
ways.  The  last  term  of  (B.ll),  which  represents  “heat  conduction ,”  corresponds  to  one 
of  many  forms  defining  the  heat  flux  of  the  electron  subsystem.  The  heat  flux  Sf  of  a 
single-component  system  is  defined  from  a  classical  relation  [13]: 

S'  =  Eci'  +  Jn  +  (P.  +  n.)  •  jjn  +  S«  (B.14) 

where  Pe  and  (le  are  the  so-called  elastic  stress  tensor  and  the  viscous  stress  tensor  of 
electrons,  respectively.  This  form  is  the  basis  of  the  hydrodynamic  model  of  electron 
transfer  in  semiconductors. 

Decomposition  of  total  energy  flux  in  (B.14)  applies  well  only  to  the  case  of  a  clas¬ 
sical  electron  gas,  which  supports  the  quadratic  relationship  between  energy  and  crystal 
momentum.  Explicit  expressions  for  all  energy  flux  components  of  degenerate  electrons 
in  the  parabolic  band  structure  can  be  found  from  the  expression  for  the  energy  current 
given  in  (B.10)  with  ae  =  0.  For  the  heat  flux,  however,  a  phenomenological  relation 
S?  cc  —  VTC  has  been  used  in  the  previous  literature  [84],  [85],  [7]. 

For  a  multi-component  system,  the  electronic  heat  flux  is  defined  differently,  and  is 
equal  to  5j?n,  since  the  total  momentum  of  the  whole  system  can  be  considered  stationary 
(the  total  mass  of  conducting  electrons  being  negligible  to  the  total  mass  of  the  whole 
system).  This  electronic  heat  flux  obviously  has  a  convective  energy  flow  component 
as  stated  in  Ref.  [2].  However,  it  is  shown  in  Sec.  B.2  that  even  the  heat  flux  Sf  for 
the  single-component  system  as  in  (B.14)  may  have  a  convective  component  of  sizable 
magnitude  which  arizes  from  particle  diffusion  in  addition  to  the  energy  flow  from  the 
temperature  gradient. 
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B.2  Heat  Flux  in  the  Hydrodynamic  Model  with  the  Relaxa¬ 
tion-Time  Approximation 

To  obtain  expressions  for  carrier  and  energy  fluxes  in  a  binary  system  consisting  of 
electrons  and  lattice,  we  have  approximated  the  collision  term  of  the  BTE  as 

(9fk/dt)c  =  -  (A  -  r)  M,  (B.15) 

so  that  the  relaxation  time  measures  the  characteristic  time  within  which  /( r,  k)  goes 
to  a  stationary  distribution,  /st,  at  a  certain  space-momentum  coordinate  (r,  k).  We 
then  have  equated  /st  with  fp(Tc),  the  nondisplaced  Fermi  distribution. 

Note  that  in  hydrodynamics  of  a  gas  of  identical  particles  which  mainly  interact  with 
themselves,  the  distribution  function  becomes  approximately  a  displaced  Maxwellian. 
This  approximation  is  valid  for  an  electron  gas  in  which  interaction  between  electrons 
and  the  crystal  lattice  is  turned  off,  since  electron-electron  scattering  will  randomize 
the  velocity  distribution  around  the  mean  velocity.  The  equation  obtained  by  inserting 
the  displaced  Maxwellian  to  /st  is  not  a  usual  form  of  relaxation-time  approximation. 
However,  the  situation  is  a  typical  example  in  hydrodynamics,  and  there  exists  an  asymp¬ 
totic  method  of  solution  to  a  general  BTE  due  to  Enskog.  According  to  this  method, 
it  has  been  shown,  up  to  the  second-order  of  Enskog’s  solution,  that  the  heat  flux  does 
not  have  a  convective  component  of  energy  flow,  and  the  first-order  expression  becomes 
simply  S®  =  -KCVTC  [87,  p.  122]. 

However,  the  real  situation  of  electrons  in  a  solid  is  a  two-  or  three-component  system 
whose  constituents  have  drastically  different  unit  masses,  so  that  the  crystal  lattice  is  vir¬ 
tually  at  rest.  When  other  forms  of  scattering  processes  (e.g.,  electron-phonon  scattering, 
ionized-impurity  scattering)  are  dominant  over  the  electron-electron  scattering,  we  are 
led  to  think  that  the  scattering  will  tend  to  randomize  the  velocity  distribution  of  each 
component  around  the  overall  mean  velocity  of  the  total  system.  We  may  suppose  that 
the  real  distribution  function  should  be  perturbed  from  the  nondisplaced  Maxwellian,  so 
that  one  can  use  fp(T)  for  />l  in  (B.15). 
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The  electron-phonon  scattering  is  indeed  dominant  over  electron-electron  scattering 
in  a  modestly  doped  semiconductor.  In  a  degenerate  semiconductor,  electron-electron 
scattering  accompanies  and  enhances  the  rate  of  other  scattering  processes  (mostly  elec¬ 
tron-phonon  scattering),  making  a  second-order  quantum-mechanical  process  [88].  The 
effect  of  high  population  makes  it  hard  to  tell  whether  the  maximum  distribution  is 
displaced  from  the  origin  k  =  o  from  the  shape  of  the  distribution  function.  Actually, 
extreme  degeneracy  (e.g.,  in  a  metal)  makes  the  following  discussion  irrelevant,  as  is 
evident  when  we  obtain  the  expression  for  the  heat  flux.  Second,  the  relaxation  time  is 
hardly  a  constant  over  different  energies,  and  the  consequent  effect  on  the  resulting  heat 
flux  component  in  the  energy  flow  in  (B.14)  comes  not  from  higher  than  first  order  in 
Tfe,  but  from  the  first  order  in  r^.  The  abnormality  in  distribution  in  k-spa.ce  caused  by 
interaction  of  electrons  with  the  crystal  lattice  gives  heat  flux,  while  nonuniformity  of 
velocity  distribution  in  real  space  gives  viscosity.  Therefore,  the  viscosity  component  of 
energy  flux  is  usually  negligible  in  electron  transport  except  in  a  region  of  high  transversal 
velocity  difference  as  in  an  inversion  layer  of  metal-oxide-semiconductor  structure  [89], 
for  instance. 

Since  the  stress  tensors  are  defined  as 

p- + ■  1  (*-£)*(*- (BI6) 

we  can  then  obtain  an  expression  for  the  heat  flux  defined  in  (B.14)  from  the  first-order 
solution  for  /k  to  the  BTE  under  the  relaxation-time  approximation.  First,  for  the  case 
that  the  relaxation  time  is  velocity-independent, 

S?  =  VTC.  (B.17) 

For  the  case  of  a  relaxation  time  Tk  a  l»k|  \  the  remaining  energy  flux  becomes  to 
the  first  order  in  Tk 

S ?  =  -MeTc  [(nr,  -  (VEC  +  TcVr,J 
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1 


(B.18) 


+  (6^2-5ii^i)  VT, 

while  the  viscosity  component  comes  from  a  term  of  order  i\tu.  Note  that  it  is  possible 
that  the  total  kinetic  energy  flow  lags  behind  the  convective  energy  flow  with  /*  obtained 
from  the  relaxation  time  oc  Note  that  when  we  specified  the  relaxation  time 

in  the  form  (B.6)  for  the  electron-lattice  interaction,  we  also  introduced  a  form  of  en¬ 
ergy  dissipation  mechanism  to  the  lattice  system.  We  then  again  separate  this  lagging 
convective  component  of  the  heat  flux  as 

s?  =  -Tc  (2 -n  -  |7> )  j,  -K'VTC,  (B.  19) 

where  Ke  =  MtTc( 6^2  —  47j.Fi)  as  in  (B.12)  with  ae  =  0.  Note  that,  for  the  case  of 
extreme  degeneracy  as  in  a  metal,  the  lagging  convective  component  of  heat  flux  again 
vanishes,  leaving  the  heat  conduction  term  only.  This  can  be  verified  by  evaluating  271  — 
§7§  with  the  formula  given  in  (B.13).  Therefore,  as  the  carrier  density  increases,  choice 
between  /m(^)  and  fM  (k  —  for  f3t  in  the  relaxation-time  approximation  becomes 
irrelevant. 

Most  previous  simulations  on  heat  flow  in  semiconductor  devices  have  been  based  on 
the  expression  for  the  electron  energy  flux  given  in  the  form  (B.14)  and  Sf  =  —KeVTc 
with  their  varying  approximations  for  the  coefficient  I<e  [84],  [90],  [85],  [86].  Though 
not  clearly  stated  in  these  works,  this  conjecture  was  based  on  the  assumption  either 
that  the  electron  gas  in  a  solid  can  be  approximated  by  a  classical  one-component  gas  (in 
which  electron-electron  scattering  represents  the  major  scattering  mechanism)  or  that  the 
first-order  relaxation-time-approximation  estimation  on  the  heat  flux  gives  a  vanishing 
coefficient.  When  we  take  care  of  the  two-component  scattering  system,  we  get  heat  flux 
component  (lagging  in  the  case  of  cc  |®k|_1  arising  from  major  scattering  processes, 
e.g.,  acoustic  phonon  scattering)  which  is  proportional  to  the  particle  flux.  We  thus 
need  to  be  aware  of  some  subtraction  of  convective  energy  flow  after  the  addition  of  the 
pressure  component. 
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APPENDIX  C 


EXPRESSIONS  FOR  CARRIER  FLUX 

There  are  several  equivalent  expressions  for  the  carrier  fluxes  and  the  energy  flux  for 
semiconductors  with  nonuniform  band  structure.  Expressions  (1.15)  and  (1.25)  are  based 
on  the  gradients  of  rje  (j ft)  with  those  of  Ec  and  Tc.  We  may  even  express  the  fluxes  in 
terms  of  the  gradients  of  Fe  and  Tc  as  shown  in  Sec.  B.2.  With  ae  =  0, 

je  =  —Me  {Fo(T)*)VFe  +  [2 Fxirje)  -  1/«*b(lfe)]  VTC} .  (C.l) 

This  form  needs  further  transformation  before  the  Scharfetter-Gummel  scheme  can  be 
applied  [24]. 

Another  type  of  expression  is  that  of  Azoff  [6]  shown  here  for  the  nonuniform  tem¬ 
perature  distribution: 

3 

-  ?ie  =  ^eV7| Tcn  +  HeTiVEc  -  -/xen7|rcVlnme.  (C.2) 

This  form  explicitly  accounts  for  the  force  term  due  to  an  effective  mass  variation  in 
space.  He  derived  this  expression  from  the  momentum  conservation  relation  with  a 
momentum-relaxation  time,  which  is  assumed  independent  of  energy.  The  suitable  form 
for  the  Scharfetter-Gummel  discretization  scheme  then  becomes 

3  e  =  -j  {liTcVn  +  n  [v  (ec  +  7  iTc)  -  \liTcV  lnme] }  .  (C.3) 

One  can  readily  construct  a  formula  similar  to  (1.27).  This  expression  can  be  interpreted 
easily  in  a  conventional  way  as  we  classify  the  flux  components:  the  diffusion  part  and 
the  drift  part.  The  drift  part  comes  not  only  from  the  conventional  force  —VEc,  but 


also  from  the  negative  gradient  of  the  electron  pressure — the  median  level  of  the  electron 
internal  energy — and  an  unusual  nonconservative  force  from  the  nonuniformity  of  the 
band  structure.  This  last  term  does  not  appear  if  we  allow  other  forms  than  Vn  as 
the  diffusion  term  as  in  (1.25).  Note  that  when  the  real  relaxation  time  is  dependent 
on  the  kinetic  energy,  we  cannot  use  the  momentum  conservation  relation  itself  for  the 
expression  of  carrier  fluxes.  Instead,  we  should  evaluate  the  moment  of  t^k  for  the  BTE. 
Then  we  will  get  different  coefficients  and  7  factors.  If  we  extend  Azoff’s  expression  for 
the  case  of  velocity-dependent  relaxation  time  t*  =  r/  |®k|  as  in  (B.6),  we  obtain 

jc  =  \  —  V  [m*M'TcTl{Th))  +  \  —  meMefo(’/e)V£c 

u  fTlg  u  TTlg 

-  |  —  m.Mere^i(*?«)Vlnme.  (C.4) 

The  choice  of  formulas  for  discretization  is  intimately  coupled  with  the  choice  of 
independent  variables  for  the  system  of  equations.  Among  the  options  are  Fe,  r]e  (or 
equivalently  Tq (rje)),  and  n.  (The  so-called  Slotboom  variables  are  not  well  suited  here.) 
Because  of  the  wide  decimal  range  of  the  carrier  densities,  the  Scharfetter-Gummel  dis¬ 
cretization  scheme  appears  to  be  a  requirement.  From  the  following  considerations,  we 
prefer  the  set  of  flux  equations  given  in  (1.25),  whose  counterpart  flux  equations  are  the 
kinetic  energy  fluxes.  First,  we  can  put  all  the  explicit  temperature  dependencies  in  less 
than  three  terms.  Second,  using  the  Onsager-symmetric  flux  equations  as  in  (1.15)— (1.20) 
makes  the  formulation  simple  and  symmetric.  The  discretization  scheme  as  introduced 
in  Sec.  1.5  is  relatively  simple  with  our  choice  of  flux  equations  (1.25)-(1.26).  Choosing 
7e  and  t/h  as  two  independent  variables  gives  considerable  flexibility  in  handling  terms 
arising  from  incomplete  ionization  of  impurities  and  from  various  recombination  rates, 
especially  the  radiative  recombination  rate  in  the  quantum-well  region. 

Note  the  difference  between  the  use  of  conventional  mobilities,  \lk  and  /Xh,  and  the 
parameters  used  here,  Me  and  M h.  The  theoretical  definition  of  mobility  is  given  in 
various  textbooks  and  is  a  dyadic  by  nature.  It  is  a  function  of  band  structure  and  the 
scattering  mechanism  through  the  scattering  time.  Use  of  any  theoretical  expression  for 


63 


these  parameters  is  usually  discouraged  in  favor  of  phenomenological  expressions  which 
depend  on  the  carrier  densities,  impurity  densities,  and  temperature.  Thus,  there  is  no 
restriction  against  using  a  new  expression  for  such  parameters  as  long  as  we  can  transform 
those  available  experimental  expressions  into  ones  that  conform  to  our  definition  by  such 
relations  as  (1.17).  In  fact,  the  new  mobility  parameters  with  Fq (v)  actually  more  closely 
represent  the  situation  in  a  semiconductor  device  in  its  degenerate  state.  Moreover,  the 
new  parameters  are  often  less  temperature  dependent  than  the  conventional  ones  [91]. 

For  the  Newton  iteration,  it  is  critical  to  have  analytic  expressions  for  the  Fermi-Dirac 
integrals  (^(r/)),  their  derivatives,  and  the  inverse  function  for  at  least  one  of  those 
Fermi-Dirac  integrals.  Availability  of  closed-form  expression  for  the  zeroth-order  Fermi- 
Dirac  integral,  To{rj)  =  ln(e”  +  1)  and  its  inverse  function  Pq1{N)  =  ln(e-/v  —  1)  deserve 
mention  with  the  special  property  j-fjiv)  =  Tj- i(v)  f°r  our  emphasis  on  integer-order 
integrals  for  flux  expressions.  Available  analytic  expressions  for  various-order  Fermi-Dirac 
integrals  were  reviewed  in  [92]  with  their  analytic  properties. 
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